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Executive  Summary 

Applied  Research  Associates,  Inc.  is  developing  a  suite  of  biomathematical  models  that  de¬ 
scribe  the  health  effects  of  nuclear  and  radiological  environments  for  scenarios  which  include 
nuclear  detonations.  In  such  scenarios,  many  radiation  injuries  combined  with  acute  injuries, 
such  as  trauma  and  thermal  burns,  are  anticipated.  Since  radiation  and  burn  have  a  sig- 
nihcant  effect  on  hematopoiesis,  we  have  developed  a  set  of  models  to  predict  the  impact 
of  these  insults  on  hematopoietic  cell  kinetics.  Modeling  these  essential  biological  processes 
will  help  improve  casualty  predictions  and  the  evaluation  of  combined  radiation  and  burn 
injury. 

In  this  study,  we  incorporate  the  effects  of  thermal  injury  into  existing  radiation  models 
of  hematopoiesis.  Models  that  describe  the  effects  of  acute  radiation  exposure  on  throm- 
bopoiesis,  granulopoiesis,  and  lymphopoiesis  were  detailed  in  a  previous  report  (Wentz  et  ah 
2014a).  Experimental  data  in  rodents  and  observational  data  in  humans  were  analyzed  to  de¬ 
termine  the  specihc  biological  processes  affected  by  burn.  Based  on  this  data,  we  hypothesize 
the  following  hematopoietic  effects  of  thermal  injury: 

•  thrombopoiesis:  decreased  platelet  lifespan,  increased  levels  of  circulating  stimulatory 
mediators  which  lead  to  an  increased  rate  of  platelet  progenitor  repopulation,  and 
increased  maturation  rate  of  immature  megakaryocytes. 

•  granulopoiesis:  increased  release  rate  of  granulocytes  from  the  bone  marrow,  de- 
margination  of  granulocytes  in  the  blood,  increased  uptake  rate  of  granulocytes  into 
the  tissues,  and  increased  rate  of  granulocyte  progenitor  repopulation. 

•  lymphopoiesis:  decreased  lymphocyte  lifespan. 

These  effects,  quantihed  as  functions  of  burn  size,  were  incorporated  into  the  radiation 
injury  models  to  create  biomathematical  models  of  radiation  combined  injury.  The  models 
accurately  describe  blood  cell  kinetics  following  thermal  injury  and  provide  predictions  of 
blood  cell  kinetics  following  combined  injury.  Compared  to  their  predictions  for  single  injury 
types  (radiation  alone  or  thermal  injury  alone),  the  models  predict  that,  in  some  cases, 
combined  injury  leads  to  lower  cell  count  nadirs  and  longer  recovery  times.  These  models 
aid  in  casualty  predictions  and  provide  insight  into  the  time-course  of  injury  and  recovery, 
medical  resource  planning,  and  disaster  preparedness. 
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Introduction 


1.1  Biomat hematical  Modeling  to  Support  DTRA  Mission 

Applied  Research  Associates,  Inc.  (ARA)  has  been  tasked  by  the  Defense  Threat  Reduction 
Agency  (DTRA)  to  support  its  mission  to  safeguard  the  United  States  against  weapons  of 
mass  destruction.  ARA  supports  this  effort  by  developing  state-of-the-art  mathematical 
models  that  predict  medical  and  performance  effects  of  radiation  and  burn  injuries,  thereby 
enhancing  our  understanding  of  the  potential  impact  of  a  nuclear  detonation.  We  are  devel¬ 
oping  mechanistic,  biologically-motivated  models  that  describe  the  link  between  insults  and 
the  probability  of  mortality. 

1.2  Modeling  Hematopoietic  Cell  Kinetics 

Radiation  and  burn  injury  lead  to  significant  alterations  in  hematopoiesis.  Cytopenias  and 
cytoses,  states  of  decreased  and  increased  blood  cell  concentrations,  respectively,  have  been 
associated  with  increased  mortality  risk  (Mirsaeidi  et  ah  2010).  Thus,  understanding  pro¬ 
genitor  and  blood  cell  kinetics  following  radiation  and  burn  can  be  useful  for  predicting 
outcomes.  We  have  developed  a  set  of  mechanistic  models  to  simulate  how  these  combined 
insults  lead  to  alterations  in  blood  cell  kinetics.  Models  that  describe  the  effects  of  acute 
radiation  exposure  on  thrombopoiesis,  granulopoiesis,  and  lymphopoiesis  were  presented  in  a 
previous  report  (Wentz  et  al.  2014a).  In  the  present  report,  we  incorporate  the  hematopoietic 
effects  of  thermal  injury  into  the  previously  developed  radiation  models  by  mathematically 
interpreting  burn-dependent  mechanisms  and  developing  new  model  parameters.  The  re¬ 
sulting  models  predict  the  effects  of  radiation  and  burn  combined  injury. 

Experimental  data  in  rodents  and  observational  data  in  humans  were  analyzed  to  deter¬ 
mine  the  specific  biological  processes  affected  by  burn.  Based  on  this  data,  we  hypothesize 
the  following  hematopoietic  effects  of  thermal  injury: 

•  thrombopoiesis:  decreased  platelet  lifespan,  increased  levels  of  circulating  stimulatory 
mediators  which  lead  to  an  increased  rate  of  platelet  progenitor  repopulation,  and 
increased  maturation  rate  of  immature  megakaryocytes  (MKs). 

•  granulopoiesis:  increased  release  rate  of  granulocytes  from  the  bone  marrow,  de- 
margination  of  granulocytes  in  the  blood,  increased  uptake  rate  of  granulocytes  into 
the  tissues,  and  increased  rate  of  granulocyte  progenitor  repopulation. 

•  lymphopoiesis:  decreased  lymphocyte  lifespan. 

These  effects,  quantihed  as  functions  of  burn  size,  were  incorporated  into  the  radiation 
injury  models  to  create  mechanistic  models  of  radiation  and  burn  combined  injury.  To 
our  knowledge,  these  are  the  first  mathematical  models  describing  the  effects  of  burn  on 
hematopoiesis.  This  report  presents  an  updated  version  of  ARA’s  thrombopoiesis  model 
that  was  previously  published  (Wentz  et  al.  2014b). 

One  of  the  most  important  advantages  of  these  biomathematical  models  is  that  they 
provide  information  that  aids  in  predicting  time-dependent  casualty  streams  and  medical 
resource  requirements  over  time.  Through  the  development  of  mechanistic  combined  injury 
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models,  we  can  gain  insight  into  the  underlying  pathophysiology  of  combined  injury,  predict 
potentially  synergistic  effects,  and  investigate  the  effects  of  countermeasures.  In  contrast, 
probit  models  do  not  permit  this  type  of  exploration. 

1.3  Implementation  of  Models 

The  models  described  in  this  report  have  been  implemented  in  Health  Effects  from  Nuclear 
and  Radiation  Environments  (HENRE),  a  user-friendly  software  platform  with  which  one 
can  (1)  dehne  an  insult,  (2)  run  models  to  predict  health  effects  of  that  insult,  and  (3)  ana¬ 
lyze  model  output.  Eventually  HENRE  will  be  integrated  with  DTRA’s  Hazard  Prediction 
and  Assessment  Capability  (HPAC),  a  modeling  platform  that  predicts  nuclear  environment 
information  such  as  radiation  dose  and  thermal  fluence  at  specific  locations.  The  environ¬ 
ment  parameters  obtained  from  HPAC  can  be  mapped  to  the  environment  parameters  that 
drive  the  combined  injury  models  described  in  this  report.  A  future  goal  will  be  to  predict 
the  effect  of  these  environments  on  heterogeneous  populations. 
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Methods 


2.1  Model  Parameterization 

To  optimize  parameters,  quantify  parameter  uncertainty,  and  assess  model  identifiability  we 
used  nonlinear  least  squares  curve  fitting  and  Markov  Chain  Monte  Carlo  (MCMC)  sam¬ 
pling.  All  ordinary  differential  equation  (ODE)  simulations,  optimization  procedures,  and 
parameter  confidence  interval  calculations  were  performed  using  the  statistical  computing 
environment  R  (vS.O.O,  R  Core  Team  2013). 

Parameter  estimation  can  be  formulated  as  a  nonlinear  least  squares  optimization  prob¬ 
lem,  and  the  structure  of  the  least  squares  objective  function  can  be  leveraged  to  create 
alternatives  to  Newton’s  method,  such  as  the  Levenberg-Marquardt  Method  (More  1978). 
The  Levenberg-Marquardt  Method,  which  can  be  viewed  as  a  combination  of  the  Gauss- 
Newton  Method  and  the  Steepest  Descent  Method,  provides  a  robust  and  efficient  way  to 
locate  optimum  parameter  values.  An  optimized  parameter  set  ^lm  was  obtained  using  the 
Levenberg-Marquardt  algorithm  implemented  in  the  modCost  and  modFit  tools  from  the 
Flexible  Modeling  Environment  (FME)  package  (Soetaert  et  ah  2010a).  As  part  of  this  pro¬ 
cess,  solution  curves  for  the  hematopoiesis  model  equations  (see  Appendix  A)  were  calculated 
using  the  “ode”  function  from  the  deSolve  package  (Soetaert  et  al.  2010b). 

ModFit  performs  constrained  fitting  of  a  model  to  data  using  input  including  (1)  a 
function  to  be  minimized,  (2)  initial  parameter  values,  (3)  the  lower  and  upper  bounds  to 
be  used  for  each  parameter,  and  (4)  the  optimization  method.  The  objective  function  /  we 
provided  to  modFit()  maps  a  given  parameter  vector  to  a  residual.  This  mapping  includes 
the  following  steps: 

1.  the  ODE  simulation  is  performed  for  each  radiation  dose  and  burn  size  combination; 

2.  the  output  of  each  ODE  simulation  and  the  observed  data  are  used  as  input  for  mod¬ 
Cost,  which  calculates  and  returns  the  residuals; 

3.  the  residuals  from  each  simulation  are  combined  and  returned. 

modFit  uses  f{^k)  and  the  Levenberg-Marquardt  Method  to  generate  the  next  approximate 
parameter  vector  These  steps  are  then  repeated  to  generate  additional  approximate 

solutions  ^fc+2,  ^fc-i-3,  •  •  •  until  a  convergence  test  is  passed.  Finally,  we  note  that  the  opti¬ 
mizations  were  performed  on  normalized  data  in  which  each  data  value  was  divided  by  a 
baseline  value.  This  allows  for  comparison  of  residuals  between  datasets  without  inadver¬ 
tently  weighting  one  dataset  more  than  another. 

Next,  using  the  parameter  set  ^lm  to  initialize  a  Markov  Chain,  the  MCMC  method  with 
the  adaptive  Metropolis  algorithm  was  used  to  estimate  conhdence  intervals  for  parameters 
(Laine  2008;  Haario  et  al.  2001).  The  MCMC  sampling  procedure  included  the  following 
steps:  (1)  given  an  initial  parameter  vector  ^j,  an  estimate  for  the  parameter  covariance 
matrix,  and  the  same  function  /  as  described  above,  the  parameter  space  was  randomly 
sampled  to  generate  a  new  parameter  vector  (2)  the  likelihood  L  of  the  new  parameter 
vector  was  determined,  and  if  L(.^j+i)  >  was  chosen;  (3)  if 

was  selected  with  probability  a,  where  a  =  P{^j+i)/P{^j)  and  P  is  the  probability 
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of  a  given  parameter  vector.  These  calculations  were  performed  in  R  using  modMCMC,  a 
function  in  the  FME  library. 

2.2  Optimization  and  Validation  Data 

To  model  the  effects  of  thermal  injury  on  blood  cell  populations,  human  and  murine  data 
were  collected  from  published  literature.  Data  on  the  time-course  of  blood  cell  concentrations 
following  thermal  injury  were  digitized  and  reviewed  (March  et  ah  2013;  Baxter  1974;  Pavic 
et  al.  2007;  Simon  et  ah  1977;  Peterson  et  ah  1983;  Vindenes  et  al.  1995;  Nijsten  et  ah  1991; 
Balch  1963;  D’Arpa  et  al.  2009;  Kagan  et  al.  1989;  Neilan  et  al.  1977).  In  some  cases,  the 
data  available  for  a  particular  cell  line  was  limited,  leading  to  broader  inclusion  criteria. 
Data  not  included  in  the  optimization  were  used  to  validate  the  resulting  models.  During 
the  optimization  procedure,  data  was  divided  by  predetermined  baseline  cell  concentrations, 
specihc  to  each  cohort  of  subjects.  The  comparison  of  data  to  the  dimensionless  form  of 
the  model  allows  for  differences  between  the  data  and  the  model  to  be  comparable  across 
studies. 

2.3  Comparisons  of  Models  with  Data 

Comparisons  of  model  simulations  with  human  data  are  challenging  due  to  the  sources  of 
variability  within  a  population.  In  the  literature,  blood  cell  kinetic  data  is  usually  provided  as 
averaged  values  across  a  group  containing  individuals  with  different  burn  sizes  and  different 
baseline  blood  cell  concentrations.  Rather  than  provide  burn  sizes  for  each  individual,  studies 
provide  the  mean,  median,  or  range  of  percent  total  body  surface  area  affected  by  burn 
(%TBSA).  In  the  models  presented,  burn  depth  is  not  currently  accounted  for.  However,  it 
likely  influences  the  hematopoietic  response,  leading  to  variability  between  patients  with  the 
same  burn  size.  For  mathematical  purposes  the  burn  size  is  denoted  as  S  in  the  equations 
and  hgures  of  this  report.  S  represents  the  burn  size  fraction  and  ranges  between  0  and  1, 
where  S  =  1  corresponds  to  100%  TBSA.  Baseline  values  for  subjects  in  human  case  studies 
are  not  available  and  must  be  assumed  based  on  distributions  in  the  general  population.  To 
understand  the  effect  of  the  variability  of  burn  size  and  blood  cell  concentration,  multiple 
model  outputs  are  compared  to  the  data  using  a  range  of  inputs. 
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3  Thrombopoiesis  Model 

Burn-dependent  parameters  were  developed  and  integrated  into  our  previously  developed 
mathematical  model  of  thrombopoiesis.  The  previous  model  developed  by  ARA  simulates 
the  effects  of  acute  radiation  exposure  (Wentz  et  ah  2014a).  Briefly,  the  existing  throm¬ 
bopoiesis  model  consists  of  three  compartments  containing  mitotic  progenitors  in  the  bone 
marrow,  post-mitotic  progenitors  in  the  bone  marrow,  and  platelets  in  circulation.  These 
compartments  are  regulated  through  feedback  mechanisms.  Acute  radiation  exposure  leads 
to  damage  and  the  eventual  death  of  a  dose-dependent  proportion  of  radiosensitive  cells.  A 
structural  diagram  of  the  thrombopoiesis  model  with  the  effects  of  burn  is  given  in  Figure  3.1. 

3.1  Modeling  the  Effect  of  Burn 

The  effects  of  burn  on  thrombopoiesis  were  incorporated  into  the  model  through  the  analysis 
of  observational  and  experimental  data  in  rodents  and  humans.  To  improve  the  model- 
selection  process,  a  model  was  hrst  developed  for  rodents  and  then  translated  to  humans. 
After  thermal  injury,  rodents  and  humans  exhibit  an  immediate  thrombocytopenia  followed 
by  recovery  and  sustained  thrombocytosis  (Kalmaz  et  al.  1991;  March  et  al.  2013).  The  likely 
mechanisms  for  the  observed  trends  are  discussed  and  incorporated  into  the  mathematical 


Cell  transitions  -  ■  Instant  radiation  damage 


Cell  transitions  affected  by  burn  . ■>  Feedback  effects 

Figure  3.1:  Thrombopoiesis  model  diagram  (adapted  from  Figure  4.1,  Wentz 
et  al.  2014a). 
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model  of  thrombopoiesis.  These  mechanisms  include  both  a  decreased  platelet  lifespan  and 
an  increased  mitotic  progenitor  cell  repopulation  rate.  The  equations  for  the  mathematical 
model  for  thrombopoiesis  with  radiation  and  burn  effects  are  given  in  Appendix  A.l. 

3.1.1  Reduction  in  Platelet  Lifespan 

Several  mechanisms  could  cause  the  immediate  thrombocytopenia  observed  following  burn. 
Various  burn  treatments  are  known  to  cause  thrombocytopenia,  but  burn-induced  throm¬ 
bocytopenia  is  also  observed  in  controlled  murine  experiments  in  which  treatment  is  not 
administered  (Kalmaz  et  ah  1991;  Fujimi  et  ah  2006;  Davis  et  ah  1955).  A  decreased  pro¬ 
duction  rate  of  platelets  could  cause  thrombocytopenia;  however,  an  increased  rate  of  platelet 
production  as  well  as  an  increased  MK  count  were  observed  in  the  early  postburn  period  in 
rodents  (Wallner  et  ah  1984;  Burenins  et  al.  1972).  This  evidence  implies  that  the  observed 
thrombocytopenia  following  burn  is  not  due  to  bone  marrow  dehciency  or  treatment. 

Experimental  and  observational  studies  have  shown  that  thermal  injury  decreases  the 
lifespan  of  circulating  platelets  in  rodents  (Burenins  et  al.  1972)  and  humans  (Simon  et  al. 
1977).  The  mechanism  of  platelet  clearance  from  the  blood  in  rodents  at  1  hour  postburn 
was  determined  to  be  extrinsic;  factors  in  the  serum  led  to  the  removal  of  platelets  from 
circulation  (Burenins  et  al.  1972).  In  this  same  study,  the  fractional  proportion  of  platelets 
at  the  burn-wound  site  increased  compared  with  controls,  while  there  was  no  fractional 
change  in  the  percent  of  platelets  in  the  spleen,  liver,  kidney,  or  lung.  Taken  together,  this 
evidence  suggests  that  the  burn-induced  thrombocytopenia  is  caused  by  the  activation  of 
platelets  and  subsequent  coagulation  cascade  at  the  burn-wound  site.  Thus,  in  our  model, 
consumption  of  platelets  is  assumed  to  be  the  sole  mechanism  of  thrombocytopenia. 

To  incorporate  this  biological  effect  into  the  mathematical  model  of  thrombopoiesis,  the 
rate  of  platelet  consumption  is  determined  as  a  function  of  burn  size  and  time  postburn. 
The  burn  insult  initially  leads  to  an  increased  rate  of  platelet  decay.  However,  as  the  time 
postburn  increases,  the  platelet  decay  rate  gradually  returns  to  the  baseline  value.  In  reality, 
burn  leads  to  tissue  damage  resulting  in  a  complex  inflammation  response  involving  platelet 
activation  (Van  Hinsbergh  2012).  Once  activated,  platelets  release  the  contents  of  granules 
into  the  serum  which,  through  a  positive  feedback  mechanism,  lead  to  the  activation  of 
more  platelets.  Rather  than  attempting  to  model  this  entire  cascade  of  events,  the  system 
is  simplihed,  and  the  rate  of  platelet  clearance  is  directly  related  to  burn  size  and  time 
postburn. 

To  determine  the  dependency  of  the  platelet  decay  rate  'ifj  on  burn  size  S  and  time 
postburn  t,  we  (1)  dehne  an  expected  platelet  lifespan  function  Tpit{S,t)  using  data  from 
the  literature,  and  then  (2)  relate  Tpit{S,t)  to  platelet  decay  rate  by  implementing  a  transit 
compartment  model  (-0  =  I / Tpit{S ,  t)) . 

To  determined  Tpit{S,  t),  we  hrst  examined  how  burn  size  S  affected  platelet  lifespan  at  a 
constant  time  postburn  t.  Data  on  platelet  lifespan  changes  in  humans  as  a  function  of  burn 
size  at  f  =  1  (i.e.,  the  survival  time  of  platelets  radio-labeled  at  f  =  1;  Simon  et  al.  1977) 
were  ht  to  linear,  exponential,  and  Hill  decay  (Hill  coefficient =1)  functions  (Figure  3.2). 
The  functions  implemented  all  contained  two  parameters;  thus,  least  squares  curve  htting 
was  used  to  determine  which  function  best  described  the  data  (see  Figure  3.2;  SSR=sum  of 
squared  residuals). 
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Figure  3.2:  Platelet  survival  time  as  a  function  of  burn  size. 


The  data  was  best  described  by  the  following  function: 

TMS,  1)  =  ^  ■  (1) 

where  V’o  =  0.106  d“^  and  Ks  =  0.112.  Note  that  V’o  is  the  rate  of  platelet  consumption 
under  healthy  conditions,  and  Ks  determines  the  burn  fraction  at  which  Tpu  is  half  the 
normal  value  {Tpit{Ks,  1)  =  l/2?/>o  =  ^pn(0,  l)/2). 

Next,  to  add  time  dependence  to  our  platelet  lifespan  model,  we  assume  the  rate  of 
recovery  does  not  depend  on  S  and  let  Ks  depend  on  time: 

t)  =  —  ■  —  ^  (2) 

Assuming  platelet  lifespan  is  minimal  immediately  following  burn  and  that  platelet  lifespan 
gradually  returns  to  normal  levels,  the  function  Ks{t)  should  have  a  minimum  at  t  =  0  and 
Ks  should  be  dehned  so  that  limt_,.oo  Ks{t)  =  oo.  We  therefore  dehne 

Ksit)  =  do  +  dit  (3) 

so 

=  (4) 

Since  Ks{l)  =  0.112,  do  =  0.112  —  di;  to  estimate  di,  we  used  data  obtained  by  Simon 
on  the  temporal  effects  of  burn  on  platelet  lifespan  (  Simon  et  al.  1977;  see  Fig.  3.3).  Five 
of  the  patients  examined  1  day  postburn  were  also  examined  3  and  5  weeks  following  burn. 


Solid  Line:  Simulated  fit 
to  Equation  5  (SSR=0.454) 

Data  (Mean  ±  standard  error): 
Simon  et  al.  1 977 


Figure  3.3:  Mean  platelet  survival  time  as  a  function  of  time  postburn. 


However,  each  of  these  patients  had  different  burn  sizes,  and  only  the  mean  platelet  lifespan 
for  the  group  was  provided  at  each  time  point.  The  mean  platelet  lifespan  for  the  5  subjects, 
each  with  a  specified  burn  size,  at  time  t  is  given  by 


1 


1  + 

^  do+dit 


1 


1  + 


Si 


(5) 


where  Si  is  the  fraction  of  surface  area  burned  for  subject  i  =  1, 2, ...,  5.  Figure  3.3  shows 
the  result  of  fitting  Equation  5  to  the  mean  platelet  survival  time  for  the  patient  group.  The 
solid  line  in  Figure  3.3  is  a  plot  of  Tpuit)  using  di  =  0.0202. 

Thus,  we  have  defined  an  expected  lifespan  function  Tpit{S,t),  given  by  Eq.  4.  Motivated 
by  the  study  by  Simon  (Simon  et  al.  1977),  in  which  platelet  survival  time  was  determined 
by  fitting  gamma  distributions  to  data  on  radio-labeled  platelet  concentrations,  we  relate 
Tpit{S,t)  to  platelet  decay  rate  by  dividing  the  platelet  compartment  into  subcompartments. 
Specifically,  in  our  model,  we  divide  the  platelet  compartment  into  m  subcompartments,  and 
the  transition  rate  between  subcompartments  is  where 


1 

Tpit{S,t)- 


This  model  structure  yields  platelet  compartment  transit  times  that  are  approximately 
gamma  distributed  (if  "0  is  constant,  the  transit  time  has  a  gamma  distribution). 


3.1.2  Increased  Mitotic  Progenitor  Repopnlation  Rate 

The  sustained  thrombocytosis  following  burn  suggests  that  burn  increases  the  rate  of  stem 
cell  repopnlation.  One  could  argue  that  this  effect  is  simply  the  system’s  homeostatic  re- 
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sponse  to  thrombocytopenia:  when  platelet  levels  drop,  negative  feedback  increases  the  stem 
cell  repopulation  rate,  and  overcompensation  leads  to  thrombocytosis.  However,  the  throm¬ 
bocytosis  would  be  transient;  that  is,  by  a  similar  argument,  platelet  concentration  would 
drop  again,  overcompensation  would  lead  to  thrombocytopenia,  and  the  cycle  would  continue 
as  the  system  oscillates  around  equilibrium  and  then  returns  to  a  steady  state.  Sustained 
thrombocytosis  is  therefore  not  consistent  with  negative  feedback  driven  by  platelet  concen¬ 
tration  per  se.  We  therefore  hypothesize  that  the  thrombocytosis  following  burn  is  due  to 
a  burn-induced  surge  of  a  stimulatory  mediator,  of  which  one  consequence  is  an  increase  in 
early  progenitor  cell  repopulation  rate. 

The  principle  mediator  of  thrombopoiesis  is  thrombopoietin  (TPO).  TPO  acts  as  a 
thrombopoietic  stimulator  while  platelets  and  MKs  provide  a  negative  feedback  mechanism 
by  receptor-mediated  uptake  and  destruction  of  TPO  (for  a  review  of  TPO  regulation  see 
Kaushansky  2005).  Interleukin  6  (IL-6)  stimulates  thrombopoiesis  in  vivo  in  mice  and  mon¬ 
keys  by  increasing  TPO  transcription  in  the  liver  (Ishibashi  et  ah  1989;  Asano  et  ah  1990; 
Kaser  et  ah  2001).  Following  a  dosing  regimen  of  IL-6  in  mice,  platelet  counts  increased 
signihcantly  by  day  4  in  a  dose-dependent  fashion  (Ishibashi  et  ah  1989).  In  monkeys, 
platelet  counts  increased  signihcantly  on  day  7  of  IL-6  treatment  (Asano  et  ah  1990).  This 
delay  between  treatment  and  effect  suggests  IL-6-induced  TPO  is  likely  acting  on  early 
thrombopoietic  progenitors. 

Increased  levels  of  IL-6  have  been  observed  within  the  hrst  few  hours  following  burn  in 
humans  (Carsin  et  ah  1997;  Drost  et  ah  1993;  Schluter  et  al.  1991)  and  remain  detectable 
through  at  least  5  weeks  postburn  (Drost  et  al.  1993).  The  increase  in  IL-6  levels  has  been 
correlated  with  burn  size  (Schluter  et  al.  1991;  Kowal-Vern  et  al.  1994),  although  one  study 
saw  no  relationship  between  burn  size  and  IL-6  levels  (Drost  et  al.  1993).  Increased  levels 
of  circulating  TPO  have  also  been  observed  in  burn  patients  both  with  and  without  sepsis 
(Lupia  et  al.  2009).  Although  elevated  TPO  levels  could  be  caused  by  a  thrombocytopenic 
state,  no  difference  in  platelet  counts  between  subject  groups  was  observed  by  Lupia  et  al. 
2009.  Lupia  et  al.  2009  quantihed  levels  of  TPO  in  burn  subjects  at  an  unspecihed  time 
following  burn.  In  burn  patients,  the  average  TPO  concentration  was  153.79  ±  33.07  pg 
mL“^,  and  in  septic  burned  patients,  it  was  374.93  ±  62.41  pg  mL“^.  These  values  are  com¬ 
pared  with  control  subject  values  of  68.90  ±  17.16  pg  mL“^.  In  burn  patients,  TPO  values 
were  increased  by  123%,  demonstrating  TPO  concentrations  can  greatly  fluctuate  following 
burn.  Collectively,  this  evidence  suggests  that  burn  leads  to  an  IL-6-dependent  increase  in 
circulating  TPO  levels,  resulting  in  an  increased  hematopoietic  stem  cell  repopulation  rate. 

To  model  this  effect  of  burn,  the  baseline  rate  of  mediator  production  G  is  multiplied  by 
the  factor  1  -|-  so  that  the  mediator  production  rate  increases  immediately  after 

burn  and  gradually  returns  to  normal  (6o  describes  the  magnitude  of  the  effect,  Oq  describes 
the  duration  of  the  effect,  and  Cq  describes  the  effect’s  dependence  on  burn  size). 

The  rate  of  change  of  mediator  concentration  is  thus 

^  =  G(1  +  -  H{xt  +  xi  +  02(4"  +  4)  +  (6) 

where  S  is  the  fraction  of  skin  surface  area  burned.  Using  a  quasi-steady-state  approxima- 


10 


(7) 


tion  for  J,  we  can  solve  Equation  6  for  I  and  express  the  repopulation  rate  B  as 


1  +  /S  {xf^  +  xf  +  02 +  xf)  +  03X3)  /(I  + 

Note  that  including  the  effect  of  burn  on  the  mediator  production  rate  has  not  altered  the 
maximum  rate  of  mitotic  progenitor  repopulation,  a. 


3.1.3  Increased  Immature  MK  Maturation  Rate 


The  two  thrombopoietic  effects  of  burn  modeled  so  far,  reduction  in  platelet  lifespan  and 
increased  repopulation  rate,  allow  one  to  describe  the  initial  thrombocytopenia  and  the  sub¬ 
sequent  sustained  thrombocytosis.  However,  without  further  modihcation  the  model  does 
not  accurately  predict  the  time  at  which  the  system  transitions  from  a  state  of  thrombocy¬ 
topenia  to  a  state  of  thrombocytosis.  The  issue  is  the  MK  maturation  time  in  the  model: 
data  show  that  platelet  levels  start  increasing  within  5  days  after  burn,  though  MKs  in  the 
model  require  more  than  5  days  to  mature.  We  therefore  hypothesize  that  burn  also  affects 
the  maturation  rate  of  MKs.  This  is  supported  by  data  that  show  that  TPO,  in  addition  to 
increasing  repopulation  rate,  also  affects  the  maturation  of  MKs  (Kaushansky  2005). 

For  these  reasons,  an  effect  of  burn  is  incorporated  into  our  model  of  the  maturation  rate 
of  immature  MKs.  Specihcally,  burn  leads  to  an  increased  transit  rate  through  the  imma¬ 
ture  MK  subcompartments.  In  the  model  we  developed  for  thrombopoiesis  following  acute 
radiation  exposure  (Wentz  et  ah  2014a),  the  transition  rate  between  immature  MK  subcom¬ 
partments  is  n6,  where  n  is  the  total  number  of  MK  subcompartments  (both  immature  and 
mature  MKs)  and 
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where  X3  is  the  platelet  concentration  at  equilibrium  (see  Eq.  25-26,  35-36  in  Wentz  et  ah 
2014a).  Equations  8  and  9  describe  a  negative  feedback  mechanism  in  which  platelet  concen¬ 
tration  regulates  immature  MK  maturation  rate.  To  describe  burn’s  effect  on  the  maturation 
rate  as  an  effect  driven  by  mediator  concentration  (as  opposed  to  platelet  concentration  per 
se),  we  multiply  6  in  Eq.  8  by  a  factor  that  is  independent  of  x^: 
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2(5o 


(10) 


Similar  to  the  parameters  that  describe  the  effect  of  burn  on  repopulation  rate,  bj  describes 
the  magnitude  of  burn’s  effect  on  MK  maturation  rate,  a/  describes  the  duration  of  the 
effect,  and  Cf  describes  the  effect’s  dependence  on  burn  size. 


3.2  Results  and  Discussion 

In  this  section,  we  present  optimization  and  validation  results  for  our  model  of  the  effect  of 
burn  on  thrombopoiesis. 
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Table  3.1:  Thrombopoiesis  model  optimization  and  validation  data. 


Reference 

Burn  Size  (S)* 

Number  of  Subjects** 

Use 

March  et  al.  2013 

0.0825  (0.015-0.15) 

84 

Optimization 

March  et  al.  2013 

0.225  (0.15-0.30) 

96 

Optimization 

March  et  al.  2013 

0.6  (0.30-0.90) 

64 

Optimization 

Baxter  1974 

0.20  (0-0.40) 

NA 

Validation 

Baxter  1974 

0.55  (0.40-0.70) 

NA 

Validation 

Baxter  1974 

0.85  (0.70-1.0) 

NA 

Validation 

Pavic  et  al.  2007 

0.05  (0-0.10) 

32 

Validation 

Pavic  et  al.  2007 

0.55  (0.10-1) 

36 

Validation 

Simon  et  al.  1977 

0.57  (0.24-0.89) 

10 

Validation 

*  Represents  the  fraction  of  surface  area  burned  (i.e.,  S'=0.50  corresponds  to  50%  of  the  surface 
area  burned).  Data  is  given  as  median  (min-max). 

**  Represents  number  of  subjects  at  beginning  of  study.  For  later  blood  cell  kinetic  data,  this 
number  is  likely  decreased  due  to  mortality  and  shorter  hospital  stays. 

NA  =  not  available 

3.2.1  Parameter  Optimization 

Values  of  6o,  oo,  cq,  &/,  a/,  and  c/  in  humans  were  determined  using  the  optimization  proce¬ 
dure  described  in  Section  2.1  and  data  on  platelet  concentrations  in  humans  following  burn 
(March  et  ah  2013).  Details  about  the  study  done  by  March  et  ah  are  provided  in  Table 
3.1.  When  optimizing  and  comparing  the  model  output  to  the  data,  the  median  burn  size 
was  used  as  input  into  the  model.  The  data  from  March  et  al.  data  was  used  for  parameter 
optimization  because  of  the  large  study  size,  temporal  resolution  of  the  platelet  data  (in 
particular  for  small  t),  long  duration  of  study,  and  detailed  information  for  three  separate 
groups  based  on  burn  size.  Because  the  initial  platelet  concentration  was  recorded  soon  after 
the  burn  incident,  and  all  three  burn  size  groups  had  average  initial  platelet  counts  within 
the  normal  platelet  concentration  range  for  humans  (140-440  x  lO^/pL;  Valentin  2002),  the 
initial  platelet  count  was  used  as  the  baseline  level  when  optimizing  (i.e.,  the  hrst  recorded 
platelet  count  was  used  as  the  model’s  equilibrium  platelet  level).  If  signihcant  decreases  in 
platelet  count  occurred  before  the  hrst  recorded  value,  the  model  would  underestimate  the 
platelet  counts. 

In  the  optimization  data,  the  group  with  the  largest  burn  size  showed  the  lowest  nadir 
while  the  group  with  the  smallest  burn  size  had  the  lowest  peak  in  the  thrombocytosis 
period.  The  goal  of  our  parameterization  was  thus  to  reproduce  these  trends  as  well  as 
approximately  match  the  observed  nadirs,  slopes  of  recovery,  and  thrombocytosis  peaks. 
Not  as  much  emphasis  was  put  on  matching  data  25  days  after  burn,  because  this  data  was 
obtained  from  subjects  with  longer  hospital  stays  with  potential  complications  resulting  in 
treatment  and  surgery. 

Table  3.2  gives  the  parameter  values  for  the  entire  human  thrombopoiesis  model.  Figure 
3.4  shows  the  model  overlaid  on  the  optimization  data.  The  depth  of  the  nadir  as  well  as 
the  subsequent  rate  of  recovery  are  accurately  predicted  by  the  thrombopoiesis  model.  In 
subjects  with  smaller  burns  (0.015  <  S  <  0.15),  the  subsequent  thrombocytosis  and  ensuing 
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0.015<S<0.15  0.15<S<0.30  0.30  <  S  <  0.90 


Time  following  burn  (days) 

Solid  Line:  Model  output  using  median  of  S,  using  first  data  point  as  baseline 
Shaded  Region:  Model  output  for  range  of  S,  using  first  data  point  as  baseline 
Data  (Mean):  Marck  et  al.  2013 

Figure  3.4:  Human  thrombopoiesis  model  compared  with  optimization  data. 


return  to  a  normal  platelet  concentration  is  also  accurately  predicted.  For  the  case  of  larger 
burns  {S  >  0.15),  the  data  from  Marck  et  al.  study  did  not  include  the  return  to  normal 
platelet  levels.  However,  in  those  cases,  it  is  possible  that  complications  leading  to  longer 
hospital  stays  led  to  an  increased  platelet  concentration.  Thus,  the  recovery  times  predicted 
by  our  model  for  S'  >  0.15  may  not  be  unrealistic. 

The  optimized  parameters  suggest  that  the  effect  of  burn  on  mediator  concentrations  is 
more  sustained  than  the  effect  on  the  immature  MK  maturation  rate  (i.e.,  Oq  <  «/).  Specif¬ 
ically,  approximately  4  days  are  required  for  the  burn  effect  on  the  mediator  concentrations 
to  be  reduced  by  half,  while  2  days  are  required  for  the  effect  on  the  MK  maturation  rate 
to  be  reduced  by  half.  Also,  the  value  of  the  parameter  bf,  which  determines  the  magnitude 
of  the  effect  of  burn  on  the  immature  MK  maturation  rate,  approached  the  upper  bound 
imposed  during  the  optimization,  50.  This  implies  the  immature  MK  maturation  time  in 
the  model  (i.e.,  the  time  spent  in  the  X21  compartment)  is  very  small  right  after  burn  injury. 
Although  biologically  this  is  unrealistic,  the  model  also  includes  a  fixed  maturation  time  in 
the  X2M  compartment  (see  Fig.  4.2  in  (Wentz  et  al.  2014a)).  This  result  suggests  that  burn 
may  be  affecting  more  than  just  the  first  half  the  of  MK  maturation  time. 

Simulated  platelet  trajectories  for  a  range  of  burn  sizes  are  shown  in  Figure  3.5.  For  burn 
sizes  ranging  from  10%  to  100%  TBSA,  as  burn  size  increases,  the  nadir  and  subsequent 
thrombocytosis  are  decreased  (see  Figure  3.5,  left  panel).  For  each  burn  size  group  there 
are  similar  peak  values  during  the  thrombocytosis  period,  suggesting  at  burn  sizes  greater 
than  15%,  the  effect  of  increasing  burn  size  is  minimal.  This  suggests  that  Cq  should  be  less 
than  1  (i.e.,  the  marginal  effect  of  burn  size  is  decreased  as  the  burn  size  increases.  A  small 
value  of  Co  causes  small  burns  to  result  in  a  large,  perhaps  unrealistic,  thrombocytosis  in 
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Table  3.2:  Biological  parameter  values  for  human  thrombopoiesis. 


Parameter 

Biological  Meaning 

Value 

a 

Maximum  rate  of  mitotic  progenitor  repopulation 

2.0  d-i 

7 

Rate  of  mitotic  progenitor  maturation 

0.26  d"^ 

^0 

Rate  of  MK  maturation  at  equilibrium 

0.06  d-i 

V’o 

Rate  of  platelet  decay  at  equilibrium 

0.11  d-i 

h 

Rate  of  damaged  cell  death 

1  d-i 

02 

Decay  rate  of  mediator  due  to  X2  cells  relative  to  Xi  cells 

0.006 

03 

Decay  rate  of  mediator  due  to  x^,  cells  relative  to  Xi  cells 

0.00018 

A 

Strength  of  feedback  from  platelets  on  MK  maturation 

3.5 

a 

Number  of  platelets  produced  per  MK 

3000 

^min 

Minimum  rate  of  MK  maturation 

0.048  d-^ 

^max 

Maximum  rate  of  MK  maturation 

0.08  d-^ 

n 

Number  of  MK  subcompartments 

10 

m 

Number  of  platelet  sub  compartments 

9 

D? 

Determines  fraction  of  damaged  xi  cells 

0.79  Gy 

rii 

Number  of  targets  per  xi  cell 

1 

no 

^2 

Determines  fraction  of  damaged  X2  cells 

2.77  Gy 

n2 

Number  of  targets  per  X2  cell 

4 

h 

Rate  of  damaged  cell  death 

1  d-^ 

bo 

Determines  max  burn  effect  on  repopulation  rate 

9.38 

ao 

Determines  duration  of  burn  effect  on  repopulation  rate 

0.17  d-^ 

Co 

Determines  how  burn  size  affects  repopulation  rate  change 

0.20 

bf 

Determines  max  burn  effect  on  MK  maturation 

49.85 

a/ 

Determines  duration  of  burn  effect  on  MK  maturation 

0.31  d"^ 

Cf 

Determines  how  burn  size  affects  MK  maturation  rate  change 

0.10 

do 

Determines  max  effect  of  burn  on  platelet  decay 

0.092 

di 

Determines  duration  of  effect  of  burn  on  platelet  decay 

0.0202  d-i 

Values  for  Bq,  Oq,  Cq,  6/,  a/,  and  Cf  were  determined  through  optimization  (see  Section  2.1). 
Values  for  do  and  di  are  discussed  in  Section  3.1.1. 

See  Wentz  et  al.  2014a  for  source  of  all  other  parameter  values. 


burn  sizes  of  less  than  1%  as  illustrated  in  Figure  3.5  (right  panel). 

We  acknowledge  that  our  parameter  optimization  was  based  on  a  retrospective  study 
(March  et  al.  2013)  and,  while  platelet  consumption  is  known  to  occur  after  burn,  other 
factors  may  influence  the  observed  thrombocytopenia.  The  initial  decline  in  platelet  concen¬ 
tration  could  be  caused  by  the  administration  of  resuscitation  fluids  or  by  drugs  known  to 
cause  thrombocytopenia,  such  as  silver  sulphadiazine,  heparin,  morphine,  and  paracetamol 
(March  et  al.  2013).  Heparin-induced  thrombocytopenia  is  known  to  occur,  however  the 
thrombocytopenia  following  burn  typically  is  much  shorter  than  that  caused  by  heparin. 
Finally,  the  thrombocytosis  period  could  be  artiflcially  caused  as  well  by  the  administration 
of  transfusions,  but  typically  transfusions  are  only  advised  when  the  nadir  is  below  50  ■  10^ 
platelets  per  /xL,  a  condition  that  was  not  observed  the  March  et  al.  study. 
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Left  Panel:  Model  output  with  S  ranging  from  0.1  to  1  (darker  lines  for  larger  S). 
Right  Panel:  Model  output  with  S  ranging  from  0.001  to  0.01  (darker  lines  for  larger  S). 

Figure  3.5:  Optimized  model  simulations  at  increasing  burn  size. 


3.2.2  Model  Validation 

Validation  of  the  thrombopoiesis  model  was  performed  by  comparing  model  ontpnts  to 
platelet  data  that  was  not  nsed  to  optimize  model  parameters.  For  validations,  we  nsed 
the  baseline  platelet  connt  mean  and  range  (250  ■  10^  and  140  ■  10^-440  ■  10^  cells  per  /xL, 
respectively)  observed  in  healthy  hnmans  (Valentin  2002).  Information  on  validation  data  is 
provided  in  Table  3.1,  and  Fignre  3.6  shows  the  validation  resnlts.  This  validation  provides 
evidence  that  the  model  is  able  to  match  trends  observed  in  platelet  concentration  following 
burn.  In  the  burn  studies  from  Baxter  1974  (Figure  3.6,  top  row),  the  platelet  concentra¬ 
tion  shows  a  faster  rebound  and  greater  thrombocytosis  than  predicted  by  the  model.  This 
recovery  could  be  caused  by  the  effects  of  treatment;  however,  the  study  does  not  provide 
information  on  treatment  protocols  for  the  subjects  included.  Also,  the  study  only  included 
patients  who  survived,  leading  to  the  removal  of  subjects  who  likely  had  lower  platelet  counts. 
Thus,  the  discrepancy  between  the  data  and  the  model  may  be  due  to  selection  bias.  The 
comparisons  shown  in  the  bottom  row  of  Figure  3.6  show  general  agreement  between  the 
data  and  simulations,  despite  the  large  burn  size  ranges  and  lack  of  error  estimates  for  the 
data. 
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Platelet  concentration  (  1 0^|iL 


a)  0  <  S  <  0.40 


c)  0.70  <  S  <  1 


b)  0.40  <  S  <  0.70 


d)0<S<0.10  e)0.10<S<1  f)  0.24  <  S  <  0.89 


Time  following  burn  (days) 

Solid  Line:  Model  output  using  median  of  S  and  250  •  as  baseline 

Shaded  Region:  Possible  model  output  for  range  of  S,  using  250  •  as  baseline 

Dotted  Lines:  Bounds  on  model  output  for  range  of  S  and  baseline  values  (140  -  440  • 

Data  (Mean):  Baxter  1974  (a,b,c);  Pavic  et  al.  2007  (d,e);  Simon  et  al.  1977  (f) 

Figure  3.6:  Human  thrombopoiesis  model  compared  with  validation  data. 
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4  Granulopoiesis  Model 

Burn-dependent  parameters  were  determined  and  integrated  into  our  previously  developed 
mathematical  model  of  granulopoiesis.  The  previous  model,  developed  by  ARA,  simulates 
the  effects  of  acute  radiation  exposure  (Wentz  et  ah  2014a).  Briefly,  the  granulopoiesis  model 
consists  of  four  compartments  containing  mitotic  progenitors  in  the  bone  marrow,  post¬ 
mitotic  progenitors  in  the  bone  marrow,  granulocytes  in  circulation,  and  granulocytes  in  the 
tissues.  These  compartments  are  regulated  through  multiple  feedback  mechanisms.  Acute 
radiation  exposure  leads  to  damage  and  the  eventual  death  of  a  dose-dependent  proportion 
of  radiosensitive  cells.  A  structural  diagram  of  the  granulopoiesis  model,  including  burn 
effects,  is  given  in  Figure  4.1. 

4.1  Modeling  the  Effect  of  Burn 

The  effects  of  burn  on  granulopoiesis  were  incorporated  into  the  model  through  the  analysis 
of  observational  and  experimental  data  in  rodents  and  humans.  To  help  guide  development 
of  the  human  model,  the  effects  of  burn  on  granulopoiesis  were  hrst  simulated  in  rodents. 
Murine  data  is  more  readily  available  and  can  be  obtained  from  controlled  environments 


Cell  transitions  -  ■  Instant  radiation  damage 


Cell  transitions  affected  by  burn  . ■>  Feedback  effects 

Figure  4.1:  Granulopoiesis  model  diagram  (adapted  from  Figure  5.1,  Wentz  et 
al.  2014a). 
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which  greatly  facilitates  the  identification  of  true  effects.  Once  developed,  the  structure  of 
the  model  was  conserved  across  species,  but  the  model  was  re-parameterized  with  human 
data. 

Following  burn  there  is  a  bimodal  granulocytosis  in  both  rodents  and  humans  (Palmer 
et  al.  2011;  Wallner  et  ah  1984;  McManus  1983;  Peterson  et  ah  1983;  Vindenes  et  ah  1995). 
Granulocyte  levels  rapidly  increase  but  then  quickly  decline  and  approach  normal  concen¬ 
trations  observed  in  healthy  humans.  There  is  a  second  granulocytosis  that,  although  not  as 
strong  as  the  first,  is  more  sustained.  The  mechanisms  behind  these  responses  and  the  corre¬ 
sponding  adjustments  to  the  mathematical  model  are  described  here.  The  equations  for  the 
mathematical  model  for  granulopoiesis  after  radiation  and  burn  are  given  in  Appendix  A. 2. 

4.1.1  Increased  Granulocyte  Demargination 

Following  burn,  in  both  rodents  and  humans  there  is  a  brief,  yet  rapid,  increase  in  periph¬ 
eral  granulocytes  in  the  blood.  Granulocytes  become  marginated  when  they  attach  to  the 
epithelial  lining  of  blood  vessels,  and  demargination  occurs  when  the  granulocytes  detach 
from  the  epithelium  and  begin  circulating  in  the  blood  once  more  (Summers  et  al.  2010). 
Epinephrine  studies  have  shown  that  thermally  injured  rats  have  a  reduced  proportion  of 
marginated  granulocytes  and  an  increased  total  number  of  peripheral  granulocytes  when 
compared  with  controls  (McManus  1983).  This  suggests  the  observed  granulocytosis  is  due, 
in  part,  to  granulocyte  demargination. 

The  compartment  for  granulocytes  in  the  blood,  X3,  represents  both  freely  circulating 
granulocytes  and  marginated  granulocytes.  Thus,  the  concentration  of  granulocytes  in  the 
blood  X3  can  be  represented  as  follows: 


X3  =  X3,c  +  X3,m  (11) 

where  xs^c  represents  granulocytes  in  the  freely  circulating  pool  and  X3,m  represents  granulo¬ 
cytes  in  the  marginated  pool.  Epinephrine  studies  have  shown  in  mice  approximately  1/4  of 
granulocytes  in  the  blood  are  in  the  freely  circulating  pool  while  3/4  are  in  the  marginated 
pool  (Johnson  et  al.  1995;  Van  Furth  et  al.  1986).  Thus,  at  equilibrium  =  \x3^rn-  In 
humans  the  ratio  of  marginated  to  freely  circulating  granulocytes  at  equilibrium  is  closer  to 
1,  and,  thus,  x^^c  =  J^3,m  (Summers  et  al.  2010). 

To  simulate  the  initial  granulocytosis  following  burn,  we  include  a  burn-dependent  effect 
on  the  proportion  of  marginated  granulocytes.  This  change,  however,  does  not  affect  the 
total  blood  granulocyte  concentration.  Moreover,  in  our  granulopoiesis  model  we  assume 
that  the  transition  rates  among  compartments  Xj,  i  =  1,2,  3, 4,  depend  only  on  the  total 
blood  granulocyte  concentration  X3  and  not  on  how  the  blood  granulocytes  are  partitioned, 
at  any  given  time  point,  into  circulating  and  marginated  pools,  x^^c  therefore  does  not  appear 
explicitly  in  our  differential  equations  (see  Appendix  A. 2).  However,  x^^c  is  the  quantity  that 
is  measured  in  experiments. 

Let  r{t)  denote  the  circulating  fraction  at  time  t,  so  that 

xsAt)  =  r{t)x3{t)  (12) 

and,  in  terms  of  normalized  model  variables. 
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X2,,c{t) 


XsA't) 

X^f 


X2.{t) 


r{t) 


where  f  is  the  circulating  fraction  at  equilibrium.  To  compare  model  output  to  normalized 
peripheral  granulocyte  concentration  data,  we  therefore  multiply  by  the  factor  r{t)/f. 
As  r{t)  increases,  a  larger  proportion  of  granulocytes  is  in  circulation,  and  x^A'^)  increases. 
For  example,  for  the  case  of  complete  demargination  (r  =  1)  in  humans  (f  =  0.5),  we 
multiply  Ts  by  2  to  estimate  the  normalized  peripheral  granulocyte  concentration. 

To  model  the  circulating  fraction  at  time  t,  we  want  to  choose  r{t)  so  that 

•  r{t)  is  maximized  immediately  following  burn, 

•  r{t)  gradually  returns  to  normal  levels, 

•  f  <  r{t)  <  1 
We  therefore  define 


r(f)=f(l  +  /,(A)e-“’-^),  (13) 

where  S  is  the  fraction  of  skin  surface  area  affected  by  burn  (0  <  S'  <  1)  and  burn  occurs 
at  f  =  0.  Note  that  the  constant  a,,  determines  the  duration  of  the  effect  and  the  function 
fr{S)  determines  the  magnitude  of  the  effect.  Also,  the  maximum  value  of  r{t)  should  be  1 
(i.e.,  r(0)  <  1),  so  we  must  have  that 

US)  <  A  1  (14) 

r 

Thus,  if  f  =  0.5,  fr{S)  <  1.  The  complete  form  of  /r(S')  is  given  in  Equation  18  and 
Appendix  A. 2.  In  summary,  the  non-dimensional  model  output  X3{t)  is  multiplied  by  1  -|- 
prior  to  comparison  to  normalized  peripheral  granulocyte  data. 


4.1.2  Increased  Rate  of  Bone  Marrow  Release 


The  other  hypothesized  mechanism  contributing  to  the  immediate  granulocytosis  following 
burn  is  an  early  bone  marrow  release  of  granulocytes  (Asko-Seljavaara  1974;  Eurenius  et  ah 
1973).  In  humans,  the  proportion  of  immature  cells  in  the  blood  increases  following  burn, 
meaning  the  transit  time  through  the  post-mitotic  compartment  is  temporarily  decreased 
(Volenec  et  ah  1979).  To  represent  this  change  in  the  model,  the  rate  of  granulocyte  release 
from  the  bone  marrow,  6,  is  given  by 


1  -|-  M{x'A  +  xA  + 


At)  =  A 


1  +  L{xf  +  xA  +  xiY 


■  (1  +  /2(S)e-“‘) 


(15) 


where  is  the  maximum  release  rate  in  the  absence  of  thermal  injury,  S  is  the  fraction  of 
skin  surface  area  affected  by  burn  (0  <  S'  <  1),  and  burn  occurs  at  f  =  0.  02  determines 
the  duration  of  the  effect,  and  /2(S')  determines  the  magnitude  of  the  effect.  Thus,  thermal 
injury  leads  to  an  increase  in  the  bone  marrow  release  rate  by  a  factor  of  1  -|-  /2(S)e““^*. 
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4.1.3  Increased  Rate  of  Entry  into  the  Tissues 

Thermal  injury  in  rats  leads  to  an  increased  granulocyte  influx  into  the  tissues,  presumably 
due  to  an  increased  rate  of  transition  from  the  serum  to  the  tissues  (Hansbrough  et  ah 
1996).  Another  potential  mechanism  for  the  increased  granulocyte  counts  in  the  tissue  is  a 
burn-induced  inhibition  in  granulocyte  apoptosis  (Chitnis  et  ah  1996).  However,  this  effect 
was  mediated  by  factors  in  the  serum  rather  than  intrinsically,  suggesting  a  slower  rate  of 
apoptosis  would  not  be  observed  in  the  tissues.  Also,  in  the  bimodal  granulocyte  response 
following  burn  there  is  a  sharp  decline  after  the  first  increase  in  granulocyte  count.  This 
steep  decline  suggests  an  increased  rate  of  transition  from  serum  to  tissues  is  occurring  rather 
than  a  decreased  apoptosis  rate.  Also,  granulocytes  could  become  demarginated  to  travel  to 
other  sites  in  the  body  for  entry  into  the  tissues. 

In  our  mathematical  model  of  granulopoiesis,  the  rate  of  entry  into  the  tissue,  k,  is 
adjusted  to  include  this  effect  as  follows: 

K{t)  =  Ko(l  +  /3(5')e"“3*)  (16) 

where  kq  is  the  rate  of  transition  in  the  absence  of  thermal  injury,  S  is  the  fraction  of  skin 
surface  affected  by  burn  (0  <  S'  <  1),  and  burn  occurs  at  t  =  0.  k  increases  immediately 
following  burn  and  then  returns  to  normal  levels  at  a  rate  determined  by  03.  The  function 
fsiS)  determines  the  magnitude  of  the  effect. 


4.1.4  Increased  Mitotic  Progenitor  Repopnlation  Rate 

Following  thermal  injury  in  rodents,  granulocyte  colony  stimulating  factor  (G-CSF)  levels 
increased  3.5  fold  1  to  2  days  postburn  (Shoup  et  al.  1998)  and  granulocyte  progenitors 
increased  82%  1  day  postburn  (Noel  et  al.  2002).  Although  this  evidence  comes  from  rodent 
studies,  we  assume  similar  trends  exist  in  humans.  Therefore,  we  hypothesize  the  granulo¬ 
cytosis  occurring  5-10  days  postburn  in  humans  is  due  to  an  increased  production  rate  of 
granulocytes.  These  biological  effects  were  implemented  in  the  model  by  altering  the  rate  of 
Xi  repopnlation,  B,  as  follows: 


m  = 


a 


1  +  +  e2xf  +  e^xf) 


+  /o(S)e-“"‘ 


(17) 


where  S  is  the  fraction  of  skin  surface  area  affected  by  burn,  burn  occurs  at  t  =  0,  Oq 
determines  the  duration  of  the  effect,  and  the  function  fo{S)  determines  the  magnitude  of 
the  effect.  Note  that  Equation  17  describes  the  repopnlation  rate  in  the  absence  of  radiation 
damage.  An  expression  for  B{t)  that  includes  the  effects  of  radiation  exposure  is  given  in 
Appendix  A. 2. 

The  alteration  to  the  repopnlation  rate  in  the  granulopoiesis  model  given  in  Equation  17 
is  different  than  the  effect  of  burn  on  the  repopnlation  rate  in  the  thrombopoiesis  model  (see 
Equation  7).  In  the  thrombopoiesis  model,  by  having  burn  act  on  the  mediator  production 
rate,  the  maximum  possible  mitotic  cell  repopnlation  rate  is  not  altered.  However,  in  the 
granulopoiesis  model,  thermal  injury  causes  an  increase  in  the  maximum  mitotic  repopnlation 
rate.  During  the  granulopoiesis  model  development,  adding  an  effect  of  burn  on  only  the 
mediator  production  rate  was  also  examined,  and  it  was  found  that  this  resulted  in  similar 
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model  outputs  as  the  effect  given  in  Equation  17.  Thus,  these  two  mechanisms  could  both 
explain  the  observed  granulocyte  kinetics  following  thermal  injury.  Both  mechanisms  are 
also  biologically  feasible.  The  maximum  mitotic  repopulation  rate  could  appear  to  increase 
due  to  burn-dependent  increase  in  cells  transitioning  into  the  granulopoietic  lineage.  Burn 
could  also  cause  a  faster  repopulation  of  mitotic  cells  already  in  the  granulopoietic  lineage. 
Both  of  these  mechanisms  would  be  mediated  by  cytokines  such  as  G-CSF.  Ultimately,  the 
chosen  mechanism  was  somewhat  arbitrary.  Thus,  during  future  model  development,  an 
effect  of  burn  on  the  mediator  production  rate  should  also  be  considered. 

4.1.5  Mapping  Burn  Size  to  Magnitude  of  Effect 

To  describe  the  effect  of  thermal  injury  on  granulopoiesis  we  have  implemented  four  modifica¬ 
tions  to  our  previously  developed  radiation  injury  model:  increased  demargination,  increased 
maximum  rate  of  bone  marrow  release,  increased  rate  of  entry  into  tissues,  and  increased 
rate  of  mitotic  progenitor  repopulation.  The  magnitudes  of  these  effects  /*  are  assumed  to 
depend  on  the  fraction  of  skin  surface  area  affected,  S.  Specifically,  we  assume  that 

MS)  =  (18) 

where  i  =  r,0,  2, 3.  These  functions  represent  the  magnitude  of  the  burn  effect  on  margina- 
tion  fr{S),  the  mitotic  progenitor  repopulation  rate  fo{S),  the  bone  marrow  release  rate 
/2(S'),  and  tissue  entry  rate  fsi^S).  br,  bo,  62,  and  63  are  constant  parameters  representing 
the  maximum  possible  effect.  A  and  k  are  also  constant  parameters  which  determine  the 
dose-response  curve. 

As  shown  in  Figure  4.2,  Equation  18  can  describe  a  variety  of  responses.  First,  if  fc  =  0 
or  A  =  0,  fi{S)  =  bi  (the  magnitude  of  effect  fi  is  independent  of  the  size  of  the  burn).  If 
k  =  1  and  A  — >■  00,  fi  is  approximately  a  linear  function  of  S.  For  larger  values  of  k,  fi  is 
more  clearly  sigmoidal.  Finally,  in  all  cases  ffO)  =  0  and  /*(!)  =  bi.  Values  for  A,  k,  and 
bi  were  determined  through  optimization.  Note  that  the  four  functions  ffS),  i  =  r,  0,2,3, 
differ  only  in  the  value  of  bi]  that  is,  we  used  the  same  “shape  parameters”  A  and  k  for  the 
four  thermal  injury  response  functions. 

4.2  Results  and  Discussion 

In  this  section  the  optimization  and  validation  results  for  the  burn  parameters  of  the  gran¬ 
ulopoiesis  model  are  presented. 

4.2.1  Parameter  Optimization 

To  determine  values  for  parameters  A,  k,  ar,  bo,  ao,  &2,  02,  &3,  and  03,  the  granulopoiesis 
model  was  optimized  to  data  on  granulocyte  concentrations  in  humans  following  burn,  br  was 
set  to  1  because  at  most  100%  of  granulocytes  in  the  blood  can  become  demarginated.  For 
optimization  and  subsequent  comparisons,  the  baseline  granulocyte  count  was  the  mean  value 
in  healthy  humans  (4.102  ■  10^  cells  per  /iL;  Valentin  2002).  The  model  was  optimized  using 
several  studies  (Table  4.1).  Due  to  the  limited  number  of  studies  that  provided  quantitative 
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Figure  4.2:  Possible  functions  defining  the  effect  of  burn  size  on  the  magnitude 
of  effects  (see  Equation  18). 


Table  4.1:  Granulopoiesis  model  optimization  and  validation  data. 


Reference 

Burn  Size  {S)* 

Group 

Number 

Use 

Peterson  et  ah  1983 

0.36  (0.30-0.50) 

Survivors 

16 

Optimization 

Peterson  et  ah  1983 

0.58  (0.43-0.75) 

Non- Survivors 

6 

Optimization 

Vindenes  et  ah  1995 

0.547  (0.21-0.95) 

27 

Optimization 

Nijsten  et  ah  1991 

0.38  (0.17-0.56) 

13 

Optimization 

Balch  1963 

0.40” 

1 

Validation 

*Represents  the  fraction  of  surface  area  burned  (i.e.,  5=0.50  corresponds  to  50%  of  the  surface  area 
burned).  Data  is  given  as  mean  (min-max). 

**No  range  provided  because  data  was  obtained  from  an  individual. 


data  on  granulocyte  concentrations  following  burn  in  humans,  four  of  the  five  studies  listed 
in  Table  4.1  were  used  to  optimize  parameters. 

Figure  4.3  shows  model  simulations  compared  with  data  used  for  optimization.  Table  4.2 
gives  the  parameter  values  for  the  entire  human  granulopoiesis  model.  In  the  granulopoiesis 
model,  we  are  able  to  simulate  the  bimodal  granulocytosis  observed  in  peripheral  granulocyte 
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a)  0.30  <  S  <  0.50 


b)  0.43  <  S  <  0.75 


Time  following  burn  (days) 


Solid  Line:  Model  output  using  mean  of  S  and  4.102  •  as  baseline 

Shaded  Region:  Model  output  for  range  of  S  using  4.102  •  lO^^L”^  as  baseline 
Data  (Mean):  Peterson  et  al.  1983  (a,b);  Vindenes  et  al.  1995  (c);  Nijsten  et  al.  1991  (d) 


Figure  4.3:  Human  granulopoiesis  model  simulations  compared  with  optimiza¬ 
tion  data. 


data  from  burn  subjects.  Based  on  the  optimized  parameter  values,  the  effects  of  burn  on 
the  marginated  granulocyte  pool  and  on  the  hematopoietic  stem  cell  repopulation  rate  have 
similar  durations  taking  approximately  1.4  d  for  the  burn-induced  change  to  be  reduced  by 
half.  The  change  in  the  bone  marrow  release  rate  and  transition  to  tissue  rate  are  more 
transient,  taking  0.5  and  0.1  days,  respectively,  to  be  reduced  by  half. 

Although  the  general  trends  in  the  granulocyte  response  are  matched,  there  are  some 
discrepancies  between  the  data  and  the  model  simulation.  The  data  from  Vindenes  et  al. 
1995  demonstrates  a  somewhat  larger  second  granulocytosis  than  the  model  predicts  (Figure 
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Table  4.2:  Biological  parameter  values  for  human  granulopoiesis. 


Parameter 

Biological  Meaning 

Value 

a 

Maximum  rate  of  mitotic  progenitor  repopulation 

2.0  d-i 

7 

Rate  of  mitotic  progenitor  maturation 

0.17  d"^ 

^0 

Determines  post-mitotic  maturation  rate 

0.34  d-^ 

m 

Determines  post-mitotic  maturation  rate 

0.27 

1 

Determines  post-mitotic  maturation  rate 

1.69 

Kq 

Normal  rate  of  transition  from  blood  to  tissue 

1.60  d-^ 

'ip 

Rate  of  granulocyte  decay  from  tissue 

1.05  d-^ 

Q2 

Decay  rate  of  mediator  due  to  X2  cells  relative  to  xi  cells 

2.24 

O3 

Decay  rate  of  mediator  due  to  X3  cells  relative  to  xi  cells 

8.45 

Oa 

Decay  rate  of  mediator  due  to  X4  cells  relative  to  Xi  cells 

8.45 

D? 

Determines  fraction  of  damaged  Xi  cells 

0.46  Gy 

^000 

Determines  ratio  of  weakly  damaged  to  damaged  Xi  cells 

2.46  Gy 

no 

Determines  fraction  of  damaged  X2  cells 

45.76  Gy 

no 

Determines  fraction  of  damaged  2:3  cells 

45.76  Gy 

no 

Determines  fraction  of  damaged  X4  cells 

45.76  Gy 

Hi 

Number  of  targets  per  cell  x*  (i=l,2,3,4) 

1 

Rate  of  damaged  cell  death 

1.0  d-i 

V 

Rate  of  weakly  damaged  cell  death  (when  t  >  Atae) 

0.08  d-^ 

T 

Determines  time  of  abortive  rise:  Atae  =  t  —  vD 

22.70  d 

V 

Determines  time  of  abortive  rise:  Atae  =  t  —  vD 

1.25  d/Gy 

f 

Fraction  of  X3  that  freely  circulates  at  equilibrium 

0.50 

A 

Determines  how  changing  burn  size  affects  the  system 

0.11 

k 

Determines  how  changing  burn  size  affects  the  system 

0.11 

br 

Determines  max  burn  effect  on  marginated  ratio 

1 

ttr 

Determines  duration  of  burn  effect  on  marginated  ratio 

0.50  d-^ 

bo 

Determines  max  burn  effect  on  repopulation  rate 

1.13  d-^ 

Oq 

Determines  duration  of  burn  effect  on  repopulation  rate 

0.49  d-^ 

&2 

Max  relative  burn  effect  on  bone  marrow  release 

20.00 

0,2 

Determines  duration  of  burn  effect  on  bone  marrow  release 

1.36  d-^ 

b3 

Max  relative  burn  effect  on  transition  to  tissue 

2.10 

0-3 

Determines  duration  of  burn  effect  on  transition  to  tissue 

7.67  d-^ 

Values  for  A,  k,  a^,  6oj  Oq;  ^2,  0-2,  ^3  and  03  were  determined  through  optimization  (see  Section  2.1). 
The  value  for  was  set  to  1  (at  most  100%  of  granulocytes  can  be  demarginated) . 

See  Wentz  et  al.  2014a  for  source  of  all  other  parameter  values. 


4.3c)  while  the  data  from  Nijsten  et  ah  1991  shows  a  smaller  granulocytosis  than  the  model 
predicts  (Figure  4.3d).  While  these  observations  contradict  each  other  to  some  extent,  they 
highlight  the  variability  observed  in  human  data  and  emphasize  the  quantitative  limitations 
of  this  modeling  effort.  In  the  Nijsten  et  al.  1991  study,  other  features,  such  as  the  hrst 
decline  in  granulocyte  concentration,  are  accurately  matched. 

One  challenge  with  the  human  granulopoiesis  optimization  is  the  lack  of  information  on 
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the  dose-dependency  of  the  rate  changes.  Based  on  the  optimized  parameter  values  for  A 
and  fc,  burn  size  has  practically  no  effect  on  the  model  output.  This  is  not  likely  the  case  in 
reality,  but  is  rather  a  consequence  of  the  available  data  used  to  optimize  the  models.  All 
the  data  were  obtained  from  groups  with  a  large  range  of  burn  sizes.  Based  on  the  data 
from  Peterson  et  al.  1983,  there  is  not  a  signihcant  difference  in  granulocyte  concentrations 
following  burn  between  survivors  (30-50%  TBSA  burned)  and  non-survivors  (43-75%  TBSA 
burned).  The  surviving  group  from  the  Peterson  et  al.  1983  study  represents  the  group 
with  the  lowest  mean  burn  size  (S'=0.36)  used  for  optimization.  The  possibility  exists  that 
a  maximal  effect  on  granulopoiesis  is  reached  at  a  certain  TBSA,  and  above  that  TBSA, 
no  greater  effect  is  observed.  However,  in  order  to  understand  the  response  at  lower  burn 
sizes,  data  from  individuals  with  burn  sizes  that  are  less  than  30%  TBSA  would  be  needed. 
Without  such  data,  our  model  predicts  that  even  very  minor  thermal  injury  causes  large 
changes  in  the  rates  discussed  in  Section  4.1.  One  way  to  manually  resolve  this  issue  would 
be  to  specify  in  the  optimization  a  burn  size  greater  than  zero  at  which  no  effect  is  expected. 

4.2.2  Model  Validation 

Validation  of  the  granulopoiesis  model  was  performed  by  comparing  model  output  to  gran¬ 
ulocyte  data  not  used  during  optimization.  Table  4.1  provides  information  on  the  valida¬ 
tion  data  used,  and  Figure  4.4  shows  the  validation  results  using  data  from  an  individual 
case  study  of  a  40%  TBSA  burn  exposure.  Here,  the  model  output  is  shown  both  using 
the  mean  granulocyte  concentration  observed  in  healthy  humans  (4.102  ■  10^  cells  per  juL; 
Valentin  2002)  as  well  as  the  range  of  granulocyte  concentrations  observed  in  healthy  humans 
(2.104  ■  10^  —  7.510  ■  10^  cells  per  /iL;  Valentin  2002).  In  this  validation,  the  model  is  able 
to  predict  the  granulocytosis  events;  however,  the  second  granulocytosis  period  in  the  data 


Solid  Line:  Model  output  using  S=0.40 
and  4.102  •  as  baseline 

Dotted  Lines:  Bounds  on  model  output  for 
range  of  S  and  baseline  values  (2.104  - 
7.510  ■  10^ piL-^} 

Data  (Individual):  Balch  1963 


Figure  4.4:  Human  granulopoiesis  model  compared  with  validation  data. 
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is  more  significant  than  the  model  predicts.  This  could  be  caused  by  a  variety  of  factors 
including  errors  in  the  model’s  structure  or  the  effect  of  factors  associated  with  human  data 
such  as  treatment.  The  patient  received  antibiotics  including  penicillin,  streptomycin,  and 
tetracyclin.  Chloromycetin  was  given  for  a  short  time  period  when  Aerobacter  aerogenes  was 
found  in  the  blood  stream  on  postburn  day  8. 

Infection  alone  can  lead  to  a  state  of  granulocytosis  (Summers  et  ah  2010).  Thus,  the 
granulocytosis  observed  in  burn  subjects  could  be,  in  part,  due  to  an  ensuing  infection.  In 
one  study  on  burn  victims  used  for  model  optimization,  the  non-survivors  all  developed 
sepsis,  while  only  2  of  the  16  survivors  developed  sepsis  (Peterson  et  ah  1983).  However, 
between  these  two  groups,  no  signihcant  difference  in  the  granulocyte  count  was  observed. 
This  suggests  that  the  effects  of  infection  on  the  granulocyte  concentration,  in  the  case  of 
burn  injury,  are  minimal. 


26 


5  Lymphopoiesis  Model 

Burn-dependent  parameters  were  developed  and  integrated  into  our  previously  developed 
mathematical  model  of  lymphopoiesis.  The  previous  model,  developed  by  ARA,  simulates 
the  effects  of  acute  radiation  exposure  (Wentz  et  ah  2014a).  Briefly,  the  lymphopoiesis 
model  consists  of  three  compartments  containing  mitotic  progenitors  in  the  bone  marrow, 
post-mitotic  progenitors  in  the  bone  marrow,  and  lymphocytes  in  circulation.  These  com¬ 
partments  are  regulated  through  various  feedback  mechanisms.  Acute  radiation  exposure 
leads  to  damage  and  the  eventual  death  of  a  dose-dependent  proportion  of  radiosensitive 
cells.  This  model  is  an  extreme  simplihcation  of  the  actual  lymphopoiesis  process.  In  reality 
lymphopoiesis  leads  to  the  production  of  many  different  types  of  lymphocytes  which  have 
different  dynamic  properties.  A  structural  diagram  of  the  lymphopoiesis  model,  including 
burn  effects,  is  given  in  Figure  5.1. 

5.1  Modeling  the  Effect  of  Burn 

The  effects  of  burn  on  lymphopoiesis  were  incorporated  into  the  model  through  the  analysis  of 
observational  and  experimental  data  in  rodents  and  humans.  Lymphocytopenia  is  commonly 


Cell  transitions  -  ■  Instant  radiation  damage 


Cell  transitions  affected  by  burn  . . >  Feedback  effects 

Figure  5.1:  Lymphopoiesis  model  diagram  (adapted  from  Figure  6.1,  Wentz  et 
al.  2014a). 
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observed  in  humans  within  the  first  few  days  following  burn  (Peterson  et  ah  1983;  Nijsten 
et  al.  1991;  Neilan  et  ah  1977;  Kagan  et  ah  1989;  D’Arpa  et  ah  2009).  The  goal  of  our  model 
development  is  to  accurately  predict  this  effect  of  burn.  The  equations  for  the  mathematical 
model  for  lymphopoiesis  after  radiation  and  burn  are  given  in  Appendix  A. 3. 

5.1.1  Increased  Blood  Lymphocyte  Clearance  Rate 

Lymphocytopenia,  or  a  net  decrease  in  the  population  of  lymphocytes,  occurs  after  burn 
injury.  However,  both  increases  and  decreases  occur  within  lymphocyte  subpopulations. 
Reductions  in  number  of  CD3+,  CD4+,  CD8’''  T  lymphocyte  levels  and  in  the  CD4+/CD8+ 
ratio  are  observed  after  thermal  injury,  while  the  number  of  CD19’''  B  cells  and  CD16’'' 
natural  killer  (NK)  cells  are  increased  (Entezami  et  al.  2010). 

The  decreased  lymphocyte  count  is  hypothesized  to  be  the  result  of  an  increase  in  lym¬ 
phocyte  activation  and  clearance  from  the  circulation.  Adhesion  and  activation  molecules  on 
lymphocytes  are  up-regulated  in  humans  following  thermal  injury  (Maldonado  et  al.  1991). 
Up-regulation  is  likely  caused  by  increased  levels  of  circulating  cytokines.  For  example,  inter¬ 
leukin  1  (IL-1),  which  is  known  to  increase  the  adhesiveness  of  lymphocytes  to  the  endothe¬ 
lium  (Cavender  et  al.  1986),  is  up-regulated  following  trauma  (Dinarello  1984).  Thus,  the 
hypothesized  mechanism  causing  lymphocytopenia,  is  a  burn-induced  increase  in  cytokine 
levels  which  in  turn  increases  the  rate  of  lymphocyte  clearance  from  the  blood.  A  second 
known  outcome  of  lymphocyte  activation  is  apoptosis  or  “activation-induced  cell  death” 
(Lebedev  et  al.  2003;  Teodorczyk-Injeyan  et  al.  1995).  In  our  model,  increased  lymphocyte 
clearance  and  lymphocyte  apoptosis  would  lead  to  an  alteration  in  the  same  parameter, 
The  effect  of  burn  on  lymphocyte  activation  is  potentially  correlated  with  burn  size.  Lym¬ 
phocyte  subsets  CD25+,  CD71+,  and  CD26’''  decreased  signihcantly  in  subjects  with  TBSA 
burn  between  70  and  90%,  while  in  subjects  with  TBSA  burned  between  30%  and  70%  the 
lymphocyte  levels  decreased  but  not  signihcantly  (Lebedev  et  al.  2003). 

Although  this  hypothesized  mechanism  can  explain  the  lymphocytopenia  following  burn, 
it  does  not  take  into  account  the  increased  numbers  of  B  cells  and  NK  cells  observed  following 
thermal  injury  (Entezami  et  al.  2010;  Kagan  et  al.  1989).  One  hypothesized  cause  for  the 
increase  in  B  cells  is  burn-induced  stress  hormones  leading  to  the  release  of  B  cells  from 
storage  sites  (Kagan  et  al.  1989).  However,  because  our  model  does  not  simulate  lymphocyte 
subpopulations,  our  goal  here  is  to  capture  the  net  effect.  Thus,  we  simplify  the  system  by 
not  including  this  stimulatory  effect  of  burn  on  lymphocyte  subpopulations. 

To  model  a  burn-induced  net  decrease  in  peripheral  lymphocytes,  a  burn-dependent 
change  to  the  blood  lymphocyte  clearance  rate  ■0  is  added.  The  new  equation  describing  the 
rate  of  lymphocyte  clearance  from  the  blood  following  a  burn  that  covers  S  ■  100%  of  the 
body’s  surface  area  is 


+  (19) 

where  S  is  the  fraction  of  surface  area  affected  by  burn,  'ipo  is  the  rate  of  lymphocyte  clearance 
from  the  blood  at  equilibrium,  and  burn  occurs  at  time  t  =  0.  represents  the  maximum 
relative  change  of  V’  which  is  assumed  to  occur  immediately  following  burn.  As  t  approaches 
inhnity,  'ip{t)  approaches  V’o  cind  the  system  returns  to  equilibrium.  The  speed  at  which 
returns  to  ipo  is  determined  by  03.  C3  determines  the  shape  of  the  burn  response  function. 


Table  5.1:  Lymphopoiesis  model  optimization  and  validation  data. 


Reference 

Burn  Size  (S)* 

Group 

Number  Use 

D’Arpa  et  ah  2009 

0.187  ±  0.048 

Non-septic 

5 

Optimization 

0.417  ±  0.09 

Survivor  septic 

11 

Optimization 

0.467  ±  0.25 

Non-survivor  septic  4 

Optimization 

Kagan  et  ah  1989 

0.302  ±  0.044” 

Survivor 

6 

Optimization 

0.623  ±  0.078” 

Non-survivor 

9 

Optimization 

Neilan  et  al.  1977 

0.41  (0.15-0.90) 

17 

Optimization 

Nijsten  et  al.  1991 

0.38  (0.17-0.56) 

13 

Validation 

Peterson  et  al.  1983 

0.36  (0.30-0.50) 

Survivor 

16 

Validation 

0.48  (0.43-0.75) 

Non-survivor 

6 

Validation 

Vindenes  et  al.  1995 

0.547  (0.21-0.95) 

27 

Validation 

*Data  is  given  as  mean  ±  SD  or  mean  (min-max) 

**Burn  size  given  for  multiple  postburn  periods.  The  value  provided  is  from  the  earliest  period. 


When  C3  equals  zero,  the  burn  size  has  no  effect.  As  C3  increases,  a  larger  burn  size  is  required 
for  a  given  effect  to  occur. 

5.2  Results  and  Discussion 

In  this  section,  the  optimization  and  validation  results  for  the  lymphopoiesis  burn  model  are 
presented. 

5.2.1  Parameter  Optimization 

To  determine  values  for  the  new  parameters  63,  03,  and  C3,  the  lymphopoiesis  model  was 
optimized  to  data  on  lymphocyte  count  in  humans  following  burn.  The  normal  mean  lym¬ 
phocyte  concentration  in  humans  was  used  as  the  baseline  value  for  the  data  (2.45  ■  10^/iL“^; 
Valentin  2002).  For  a  summary  of  studies  used  in  the  optimization,  see  Table  5.1. 

Figure  5.2  compares  the  model  simulation  with  data  used  for  the  optimization,  and  Table 
5.2  gives  the  parameter  values  for  the  model.  The  optimized  value  of  63  is  120.55,  which 
implies  that  a  100%  TBSA  burn  would  cause  '0  to  equal  1.22  d“^  immediately  following 
burn.  This  leads  to  an  expected  lymphocyte  transit  time  in  the  circulation  of  0.82  days 
compared  with  100  days  at  equilibrium.  The  optimized  value  of  03  is  0.86  d“^,  meaning  it 
takes  0.81  days  for  the  change  in  the  lymphocyte  clearance  rate  to  decrease  by  half. 

Given  available  data  and  the  structural  constraints  of  the  radiation-effects  lymphopoiesis 
model,  an  effect  of  burn  was  added  to  the  blood  clearance  rate  of  lymphocytes.  In  some  case 
studies  our  optimization  results  match  the  data  fairly  well: 

•  Non-septic  and  survivor  septic  groups  from  D’Arpa  et  ah  2009  (Figure  5.2ab) 

•  Non-survivor  group  from  Kagan  et  ah  1989  (Figure  5.2e) 

•  Data  from  Neilan  et  ah  1977  (Figure  5.2f) 
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a)  S=0.1 87  +/-  0.048  b)  S=0.41 7  +/-  0.09 


c)  S=0.467  +/-  0.25 


o 

o 


d)  S=0.302  +/-  0.044  e)  S=0.623  +/-  0.078 


f)0.15<S<0.90 


Time  following  burn  (days) 

Solid  Line:  Model  output  using  mean  of  S  and  2.45  •  as  baseline 

Shaded  Region:  Model  output  for  range  of  S  and  baseline  values  (1.40  -  4.20  • 

Data  (Mean  ±  SD):  D’Arpa  et  al.  2009  (a,b,c);  Kagan  et  al.  1 989  (d,e);  Neilan  et  al.  1 977  (f) 


Figure  5.2:  Human  lymphopoiesis  model  compared  with  optimization  data. 


In  other  cases  the  model  does  not  capture  the  observed  trends: 

•  Non-survivor  septic  group  from  D’Arpa  et  al.  2009  (Figure  5.2c) 

•  Survivor  group  from  Kagan  et  al.  1989  (Figure  5. 2d) 

In  the  lymphocyte  kinetic  survivor  data  from  Kagan  et  al.  1989,  the  recovery  rate  is  faster 
than  that  predicted  by  the  model  (Figure  5. 2d).  This  may  be  caused  by  selection  bias 
in  separating  survivors  from  non-survivors  with  the  same  burn  sizes.  For  the  lymphocyte 
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Table  5.2:  Biological  parameter  values  for  human  lymphopoiesis. 


Parameter 

Biological  Meaning 

Value 

a 

Maximum  rate  of  mitotic  progenitor  repopulation 

2.0  d-i 

7 

Rate  of  mitotic  progenitor  maturation 

0.06  d"^ 

Rate  of  post-mitotic  progenitor  maturation 

0.001  d-i 

V’o 

Normal  rate  of  lymphocyte  clearance  from  the  blood 

0.01  d-i 

Qx 

Mediator  decay  due  to  Xi  cells  relative  to  decay  due  to  Xi 

cells 

1 

02 

Mediator  decay  due  to  X2  cells  relative  to  decay  due  to  Xi 

cells 

2.67 

03 

Mediator  decay  due  to  X3  cells  relative  to  decay  due  to  Xi 

cells 

2.67 

D? 

Determines  fraction  of  damaged  xi  cells 

3.98  Gy 

ni 

Number  of  radiation  targets  per  cell  xi 

1 

no 

^2 

Determines  fraction  of  damaged  X2  cells 

3.98  Gy 

n2 

Number  of  radiation  targets  per  cell  X2 

1 

no 

-^3 

Determines  fraction  of  damaged  0:3  cells 

6.22  Gy 

ns 

Number  of  radiation  targets  per  cell  0:3 

0.40 

Rate  of  damaged  cell  death 

1  d-i 

h 

Determines  maximum  effect  of  burn  on  blood  clearance  rate 

120.55 

Os 

Determines  duration  of  effect  of  burn  on  blood  clearance  rate 

0.86  d-^ 

C3 

Determines  how  burn  size  affects  rate  changes 

0.84 

Values  for  63,  03,  and  C3  were  determined  through  optimization  (see  Section  2.1). 
See  Wentz  et  al.  2014a  for  source  of  all  other  parameter  values. 


kinetic  data  from  non-survivors  who  developed  sepsis  (Figure  5.2c;  D’Arpa  et  al.  2009),  the 
model  predicts  a  less  substantial  nadir  than  that  given  by  the  data.  The  development  of  sepsis 
leads  to  an  increase  in  lymphocyte  apoptosis  which  ultimately  leads  to  lymphocyte  depletion 
(Hotchkiss  et  al.  2001).  Thus,  observing  a  lower  nadir  in  the  data  is  not  unreasonable.  The 
model  optimization  would  benefit  from  the  collection  of  more  data  on  lymphocyte  counts 
following  burn  so  that  septic  patients  are  removed  from  the  optimization,  and  ideally  subjects 
would  be  grouped  according  to  burn  size  rather  than  survival. 

5.2.2  Model  Validation 

Validation  of  the  lymphopoiesis  model  was  performed  by  comparing  model  outputs  to  lym¬ 
phocyte  data  that  was  not  used  in  the  optimization.  For  these  comparisons,  the  mean 
and  range  of  the  lymphocyte  concentration  observed  in  healthy  humans  was  used  (2.45; 
1.40  —  4.20  ■  10^/rL“^;  Valentin  2002).  Table  5.1  provides  a  summary  of  the  data  used  in 
validation,  and  Figure  5.3  shows  the  results  of  these  comparisons.  These  data  were  reserved 
for  validation  as  opposed  to  optimization  due  to  the  lack  of  error  measurements  on  the  lym¬ 
phocyte  concentration  values.  Without  an  estimate  of  the  error,  it  is  difficult  to  assess  the 
predictive  quality  of  our  model.  However,  it  is  possible  to  assess  the  model’s  ability  to  match 
general  trends.  In  all  four  studies  used  for  validation,  there  is  a  decrease  in  lymphocyte  val¬ 
ues  below  the  mean  concentration  observed  in  healthy  human  populations,  which  the  model 
reliably  predicts.  However,  in  some  cases,  the  model  does  not  accurately  predict  the  time 
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a)  0.30  <  S  <  0.50  b)  0.43  <  S  <  0.75 


Time  following  burn  (days) 

Solid  Line:  Model  output  using  mean  of  S  and  2.45  •  as  baseline 

Shaded  Region:  Model  output  for  range  of  S  using  2.45  •  as  baseline 

Dotted  Lines:  Bounds  on  model  output  for  range  of  S  and  baseline  values  (1.40  -  4.20  • 

Data  (Mean):  Peterson  et  al.  1983  (a,b);  Nijsten  et  al.  1991  (c);  Vindenes  et  al.  1995  (d) 

Figure  5.3:  Human  lymphopoiesis  model  compared  with  validation  data. 


to  recovery.  Peterson  et  al.  1983  collected  lymphocyte  kinetic  data  from  survivors  and  non¬ 
survivors  (Figure  5.3a  and  5.3b,  respectively).  Our  model  predicts  a  faster  recovery  for  the 
non-survivors  than  that  seen  in  the  data  and  a  slower  recovery  for  the  survivors  than  that 
seen  in  the  data.  This  discrepancy  may  be  again  due  to  the  grouping  of  patients  into  survivor 
and  non-survivor  groups  since  lymphocyte  recovery  seems  to  be  correlated  with  survivability. 
The  models  developed  here  are  intended  to  present  the  average  response  in  humans  after 
thermal  injury,  and  therefore,  it  may  not  be  possible  to  match  outputs  to  stratified  patient 
data.  Furthermore,  sepsis  could  be  leading  to  a  delayed  recovery  in  the  non-survivor  group. 
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6  Combined  Injury  Model  Analysis 

Our  aim  in  the  current  work  is  to  develop  computational  models  that  predict  hematopoietic 
responses  of  combined  radiation  and  burn  injury  given  the  limited  quantitative  clinical  data 
available  of  combined  injury  in  humans.  To  predict  blood  and  progenitor  cell  dynamics 
following  combined  radiation  and  burn  injury,  the  effects  of  burn  were  integrated  into  pre¬ 
viously  developed  hematopoietic  models  (Wentz  et  ah  2014a).  These  single  injury  models 
were  validated  using  data  from  human  patients  who  were  accidentally  irradiated  or  thermally 
injured.  Due  to  a  lack  of  data  from  patients  with  both  radiation  and  burn  injury,  validation 
of  the  combined  injury  model  using  human  data  is  not  possible.  However,  to  evaluate  the 
credibility  of  the  combined  injury  predictions,  the  structure  of  the  combined  injury  model  is 
validated  by  comparing  murine  data  to  murine  versions  of  the  models.  In  generating  murine 
models,  the  model  structure  was  conserved  across  species  but  was  parameterized  using  only 
rodent  data.  Details  of  this  parameterization  for  thrombopoiesis  and  granulopoiesis  are  pro¬ 
vided  in  Appendix  B.  A  murine  model  simulating  the  combined  effects  of  burn  and  radiation 
on  lymphopoiesis  does  not  currently  exist,  and,  therefore,  we  were  unable  to  perform  the 
comparisons  between  the  murine  combined  injury  data  and  the  murine  model  output.  In¬ 
stead,  for  lymphopoiesis  the  murine  combined  injury  data  is  compared  to  predictions  given 
by  the  human  model. 

Blood  cell  concentrations  following  radiation  only,  burn  only,  and  combined  injury  were 
obtained  from  studies  in  mice  (Palmer  et  ah  2011;  Boudagov  et  ah  2006;  Tajima  et  ah  2013; 
Kiang  et  ah  2014),  rats  (Davis  et  ah  1955),  and  dogs  (Reid  et  ah  1955).  Since  no  canine 
models  exist,  the  canine  data  was  compared  to  human  model  outputs.  The  animal  data 
were  normalized  using  data  from  sham  controls  when  available.  Otherwise,  the  data  were 
normalized  using  the  control  blood  cell  values  prior  to  injury. 

6.1  Platelets 

Figure  6.1,  6.2,  6.3,  and  6.4  compare  radiation,  burn,  and  combined  injury  platelet  data 
in  rodents  to  murine  thrombopoiesis  model  predictions.  In  the  data  from  rat  experiments, 
shown  in  Figure  6.1  and  6.2,  the  combined  injury  and  radiation  alone  data  tend  to  follow 
similar  trajectories.  The  data  does  suggest  that  following  combined  injury  a  faster  recovery 
occurs  than  following  only  irradiation.  This  trend  towards  a  faster  recovery  is  predicted 
by  the  murine  model.  In  the  study  done  by  Boudagov  et  ah  2006  in  mice  (Figure  6.3), 
relative  comparisons  between  individual  insults  and  combined  injury  effects  are  challenging 
because  radiation  and  burn  only  data  are  not  available  past  three  days.  During  this  hrst 
3-day  interval,  the  combined  injury  response  appears  worse  than  both  the  radiation  and  burn 
only  response.  Again,  the  early  synergistic  effect  of  burn  and  radiation  is  predicted  by  the 
model.  An  exacerbated  early  response  is  a  reasonable  prediction  due  to  radiation-induced  cell 
death  and  increased  platelet  consumption  caused  by  the  reaction  to  burn.  However,  a  faster 
recovery  is  also  predictable  due  to  the  stimulatory  effect  of  burn  on  the  early  progenitors. 

Kiang  et  ah  2014  quantihed  platelet  cell  counts  at  30  days  post  9.5  Gy  and  15%  TBSA 
burn  combined  injury  in  mice  (Figure  6.4).  The  model  is  able  to  predict  the  relative  differ¬ 
ence  in  platelet  concentration  between  injury  types.  However,  compared  with  sham  values 
platelet  counts  are  decreased  following  any  combinations  of  insults  (Kiang  et  ah  2014  spec- 
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25%  TBSA 


0.96  Gy 


Cl:  25%  TBSA,  0.96  Gy 


Figure  6.1:  Platelet  concentration  following  low  radiation  dose  combined  injury 
in  rats  (top;  Davis  et  al.  1955)  compared  with  murine  model  simulations  (bot¬ 
tom). 


ifies  that  only  radiation  alone  and  combined  injury  led  to  a  significant  decrease),  while  the 
model  predicts  platelet  counts  are  close  to  normal  for  all  three  injury  combinations.  In  de¬ 
veloping  the  radiation  only  model  for  mice,  the  highest  radiation  dose  used  for  optimization 
or  validation  was  7.68  Gy  (see  Figure  B.l  and  B.2).  Thus,  the  model  may  not  currently  be 
able  to  capture  the  platelet  response  following  9.5  Gy.  In  the  burn  only  optimization  data. 
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Cl:  15%TBSA,  4.8  Gy 


15%TBSA 


4.8  Gy 


Time  following  insults  (days) 

Figure  6.2:  Platelet  concentration  following  high  radiation  dose  combined  in¬ 
jury  in  rats  (top;  Davis  et  al.  1955)  compared  with  murine  model  simulations 
(bottom). 


platelet  concentrations  were  either  normal  or  higher  than  normal  at  close  to  30  days  postbnrn 
(Fignre  B.3).  Thus,  although  Kiang  et  al.  2014  saw  a  trend  towards  a  decreased  platelet 
concentration  following  only  burn,  this  is  not  a  consistent  observation  across  all  combinations 
of  insult  ranges.  The  cause  of  the  variability  in  platelet  response  between  experiments  may 
be  due  to  differences  in  burn  depth.  This  is  not  currently  accounted  for  in  the  models. 
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Platelet  concentration  (normalized) 


10%TBSA 


7  Gy 


Cl:  1 0%  TBSA,  7  Gy 


Time  following  insults  (days) 


Figure  6.3:  Platelet  concentration  following  combined  injury  in  mice  (top; 
Boudagov  et  al.  2006)  compared  with  murine  model  simulations  (bottom). 
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Platelet  concentration  (normalized) 


1.25 

1.00- 

0.75- 

0.50- 

0.25 

0.00-1 


1 5%  TBS  A  9.5  Gy  Cl:  1 5%  TBS  A 

9.5  Gy 


Figure  6.4:  Concentration  of  platelets  30  days  following  combined  injury.  Data 
in  mice  (left;  Kiang  et  al.  2014)  compared  with  model  predictions  (right). 
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6.2  Granulocytes 

Comparisons  of  granulopoiesis  rodent  data  with  murine  model  simulations  are  shown  in 
Figure  6.5  and  6.6.  Although  there  are  deviations  between  the  data  and  the  model  predic¬ 
tions,  key  trends  are  matched  (Figure  6.5).  Exposure  to  both  radiation  and  burn  leads  to 
an  ameliorated  response  compared  with  exposure  to  only  radiation.  Specifically,  following 


15%  TBSA  (Palmer  2011) 

— ♦— 

5  Gy  (Palmer  2011) 

Cl:  1 5%  TBSA,  5  Gy  (Palmer  201 1 ) 

Cl:  1 0%  TBSA,  7  Gy  (Boudagov  2006) 
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Figure  6.5:  Granulocyte  concentration  following  combined  injury  in  mice  (top) 
compared  with  murine  model  simulations  (bottom). 
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Figure  6.6:  Concentration  of  granulocytes  30  days  following  combined  injury. 
Data  in  mice  (left;  Kiang  et  al.  2014)  compared  with  model  predictions  (right). 


combined  injury  there  is  a  less  significant  nadir  and  recovery  occurs  earlier  than  following 
only  radiation.  These  effects  are  predicted  by  the  murine  model.  Figure  6.5  also  shows  a 
comparison  of  two  combined  injury  simulations  with  different  insults.  The  model  and  the 
data  both  show  that  10%  TBSA  burn  and  7  Gy  irradiation  leads  to  a  worse  response  than 
15%  TBSA  burn  and  5  Gy  irradiation. 

In  a  second  comparison  granulocyte  concentrations  at  30  days  post  insult  are  compared 
with  model  predictions  (Figure  6.6).  The  radiation  and  combined  injury  model  outputs 
match  the  data;  however,  following  only  burn  granulocyte  levels  are  close  to  normal  while 
the  model  predicts  that  they  are  elevated.  Other  data  obtained  from  mice  with  15%  TBSA 
burn  shows  the  reverse  trend  in  which  the  model  is  underpredicting  the  platelet  response  at  30 
days  postburn  (Figure  6.5).  This  demonstrates  that  deviations  from  the  model  and  the  burn 
only  data  may  be  due  to  experimental  variability  (e.g.,  differences  in  burn  depth/severity). 

Additional  combined  injury  data  for  granulopoiesis  are  shown  in  Figure  6.7.  This  second 
set  of  data  comes  from  a  ffow  cytometry  analysis  of  neutrophils  and  other  blood  cells  (Tajima 
et  al.  2013).  Information  on  the  total  cell  population,  as  well  as  the  relative  proportion  of 
each  cell  type  in  the  population,  was  given.  However,  in  some  cases  the  total  percentage  was 
greater  than  100%,  suggesting  the  absolute  values  obtained  using  this  data  may  be  inaccu¬ 
rate.  Thus,  when  examining  this  data  we  are  mainly  interested  in  the  relative  differences 
between  combined  injury,  burn  only,  and  radiation  only  exposure.  Based  on  this  data  it  is 
again  clear  that,  as  predicted  by  the  murine  model,  adding  burn  to  a  radiation  injury  leads 
to  a  less  severe  granulocyte  blood  cell  concentration  nadir. 
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Figure  6.7:  Granulocyte  concentration  following  25%  TBS  A  burn  and  irradiation 
in  mice  (Tajima  et  al.  2013). 


6.3  Lymphocytes 

Figure  6.8  and  6.9  show  comparisons  of  the  human  lymphopoiesis  model  with  murine  data. 
The  lymphopoiesis  model  predicts  that  combined  injury  increases  the  lymphocyte  depletion 
when  compared  with  only  radiation  or  burn  injury  (Figure  6.8).  This  effect  is  suggested  by 
the  murine  data  at  later  time  points  (i.e.,  14  and  28  days  after  the  insult).  Furthermore, 
at  30  days  post  exposure,  the  model  is  able  to  match  experimental  predictions  (Figure  6.9). 
However,  at  earlier  time  points  (6  h,  24  h,  48  h),  radiation  resulted  in  a  more  significant 
nadir  than  combined  injury  (Figure  6.8).  Following  only  burn,  the  data  shows  an  initial  brief 
increase  in  lymphocyte  concentration,  which  the  model  does  not  predict.  The  lymphopoiesis 
model  is  an  extremely  simplified  version  of  the  biological  reality.  Lymphocytes  represent 
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Figure  6.8:  Lymphocyte  concentration  following  combined  injury  in  mice  (top; 
Palmer  et  al.  2011)  compared  with  human  model  simulations  (bottom). 


a  heterogeneous  population  of  cells  (e.g.,  B  cells,  T  cells,  NK  cells)  which  have  different 
dynamic  properties.  In  fact,  some  lymphocytes  are  known  to  proliferate  (Macallan  et  al. 
2005),  and  this  dynamic  is  not  currently  in  the  model  but  may  explain  the  slight  increase 
observed  in  lymphocyte  concentration  following  burn.  This  increase  may,  in  turn,  contribute 
to  the  less  severe  nadir  observed  following  combined  injury.  Therefore,  a  more  detailed  model 
of  lymphopoiesis  that  can  capture  these  effects  may  be  warranted  in  the  future. 
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Figure  6.9:  Concentration  of  lymphocytes  30  days  following  combined  injury. 
Data  in  mice  (left;  Kiang  et  al.  2014)  compared  with  model  predictions  (right). 


Figure  6.10  shows  additional  lymphocyte  combined  injury  data.  This  second  set  of  data 
comes  from  the  same  study  as  the  data  presented  in  Figure  6.7  (Tajima  et  al.  2013).  To 
quantify  the  lymphocyte  population,  we  pooled  flow  cytometry  data  for  CD4  T  cells,  CDS 
T  cells,  NK  cells,  and  B  cells.  Due  to  the  previously  discussed  limitations  of  the  study,  when 
examining  this  data  we  are  mainly  interested  in  the  relative  differences  between  combined 
injury,  burn  only,  and  radiation  only  exposure.  At  low  doses  (0.5  and  1  Gy),  there  is  no 
clear  trend  between  the  radiation  only,  burn  only,  and  Cl  data.  At  2  Gy  combined  injury 
leads  to  a  larger  decrease  in  lymphocyte  counts.  At  higher  doses  (>  4  Gy)  the  radiation 
only  and  combined  injury  response  have  similar  kinetic  profiles.  Thus,  it  appears  that  for 
higher  doses  and  a  TBSA  burn  of  25%,  radiation  is  the  dominant  insult.  Based  on  the  data 
shown  in  Figure  6.8  and  6.9,  it  appears  that  burn  size  has  an  effect  on  the  early  lymphocyte 
dynamics.  Specifically,  the  decrease  in  lymphocyte  following  combined  injury  is  less  when 
the  burn  size  is  smaller.  This  effect  of  burn  is  not  fully  understood  and  may  require  a  more 
complex  model  of  lymphopoiesis  to  accurately  simulate  the  lymphocyte  response. 
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6.10:  Lymphocyte  concentration  following  25%  TBS  A  burn  and  irradia- 
mice  (Tajima  et  al.  2013). 
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CI:20%TBSA,  0.96  Gy 


20%  TBSA 


0.96  Gy 


Figure  6.11:  White  blood  cell  dynamics  following  combined  injury  in  dogs  (top; 
Reid  et  al.  1955)  compared  with  human  model  simulations  (bottom). 


6.4  White  Blood  Cells 

Figure  6.11  shows  a  comparison  of  the  human  model’s  prediction  of  white  blood  cell  counts 
with  canine  data.  For  this  comparison  the  granulocyte  and  lymphocyte  human  model  pre¬ 
dictions  were  added  together  using  4.102  and  2.45  ■  10^/iL“^  as  the  baseline  granulocyte  and 
lymphocyte  concentrations,  respectively.  The  model  successfully  predicts  that  combined 
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injury  leads  to  an  early  ameliorated  response  in  the  white  blood  cell  concentration. 

6.5  Combined  Injury  Analysis  Summary 

Based  on  the  combined  injury  predictions  of  the  three  different  blood  cell  lines,  the  effects 
of  combined  injury  are  not  simply  exacerbated  from  what  is  observed  with  single  injuries. 
In  some  instances  combined  injury  leads  to  a  more  severe  nadir  and/or  a  delayed  recovery. 
On  the  other  hand,  the  models  suggests  that  combined  injury  could  accelerate  recovery 
relative  to  radiation  injury  alone.  The  mathematical  models  help  to  delineate  the  potential 
interactions  or  contradictory  effects  of  each  injury  when  observed  in  combination.  To  fully 
understand  the  impact  of  radiation  and  burn  combined  injury,  the  next  step  will  be  to 
correlate  features  of  the  curves  (e.g.,  nadir  and  duration  of  cytopenia)  with  risk  of  infection, 
hemorrhage,  and/or  mortality. 
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7  Future  Work 


Hematopoietic  models  describing  the  effects  of  combined  radiation  and  burn  injury  on  blood 
cell  kinetics  have  been  developed.  Additional  steps  will  be  taken  so  that  the  models  have 
direct  operational  use  for  nuclear  disaster  scenario  analysis  and  preparedness  planning. 

7.1  Validating  Combined  Injury  Model 

As  mentioned  previously,  validation  of  the  human  combined  injury  model  with  data  on  vic¬ 
tims  who  received  both  radiation  and  thermal  injury  is  ultimately  needed.  This  is  a  necessary 
step  in  verifying  that  we  are  correctly  accounting  for  the  combined  effects  of  radiation  and 
burn.  Literature  on  radiation  incidents  will  be  re-examined  in  an  effort  to  identify  human 
case  studies  involving  thermal  burns  and  radiation  exposure  that  have  blood  cell  kinetics 
data.  Some  patients  from  Nagasaki  and  Hiroshima  were  known  to  have  suffered  combined 
injury.  However,  detailed  information  on  blood  cells  kinetics  from  this  population  has  not  yet 
been  identified.  Another  mechanism  for  adding  credibility  to  our  model  predictions  would 
be  to  test  these  predictions  in  animal  experiments.  This  would  involve  a  collaborative  efforts 
with  experimentalists  actively  involved  in  combined  injury  research  and  will  be  dependent 
on  the  flexibility  of  their  ongoing  research. 

7.2  Correlating  Model  Outputs  with  Mortality 

Efforts  are  underway  to  correlate  model  outputs  (i.e.,  blood  cell  population  nadirs,  the 
duration  of  cytopenias)  with  the  probability  of  and  time  to  mortality.  If  needed,  a  weighted 
combination  of  model  outputs  may  be  used  to  assist  in  this  prediction.  Again,  correlation  of 
mortality  risk  with  clinical  endpoints  as  provided  by  the  models  can  be  guided  and  verified 
to  some  extent  by  experimental  animal  data. 

Although  we  plan  to  link  physiological  outputs  from  the  models  presented  here  to  mortal¬ 
ity  predictions,  other  hematopoietic  properties  not  predicted  by  the  models  may  contribute 
to  mortality.  Potential  future  modeling  endeavors  that  may  assist  in  mortality  predictions  are 
discussed  below.  However,  these  efforts  would  only  be  considered  if  a  satisfactory  mortality 
prediction  model  cannot  be  developed  using  the  models  as  they  stand  now. 

Thermal  injury  leads  to  a  variety  of  changes  to  the  granulopoietic  system  that  are  not 
predicted  by  our  model  as  it  currently  stands.  Studies  have  suggested  that,  following  burn, 
there  is  a  shift  of  granulocyte-monocyte  progenitors  towards  the  monocyte  lineage  (Santan- 
gelo  et  ah  2001;  Noel  et  ah  2007;  Silva  et  ah  2002).  Monocyte  concentrations  have  been 
correlated  with  survival  following  thermal  injury  (Peterson  et  ah  1983).  Thus,  monocytes 
may  prove  to  be  a  predictor  of  mortality  risk.  Future  work  may  explore  the  feasibility  of 
developing  a  model  that  can  predict  monocyte  dynamics  after  burn  and  radiation  if  adequate 
data  exists. 

Also,  neutrophil  function  has  been  shown  to  decrease  following  burn  (Arturson  1985). 
Although  this  does  not  affect  the  kinetics,  it  may  play  a  large  role  in  determining  an  indi¬ 
vidual’s  susceptibility  to  infection  following  burn  and  subsequent  risk  of  mortality.  Thus, 
future  work  might  involve  developing  a  model  to  quantify  granulocyte  dysfunction  following 
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burn.  This  could  be  done  by  creating  an  injured  granulocyte  compartment  and  determin¬ 
ing  the  proportion  of  granulocytes  that  are  injured.  Another  approach  might  be  to  allow 
granulocytes  to  have  a  range  of  injury  levels  that  follows  a  distribution  scheme.  In  turn, 
the  distribution  of  injury  levels  could  be  correlated  with  the  likelihood  of  infection.  These 
alternate  effects  of  thermal  injury  provide  other  areas  of  exploration  which  may  help  predict 
survivability  following  combined  injury. 

As  mentioned  previously,  lymphocyte  subsets  show  a  variety  of  responses  after  thermal 
injury.  However,  our  current  lymphopoiesis  model  only  predicts  the  net  effect  on  these 
subsets.  The  ratios  of  different  lymphocyte  subsets  (i.e.,  the  number  of  T  cells  compared 
to  the  number  of  B  cells)  may  be  of  importance  for  understanding  the  capabilities  of  the 
immune  system  and  whether  it  is  fully  functional.  Thus,  potential  model  developments  may 
involve  dividing  the  peripheral  lymphocyte  populations  into  various  sub-populations  so  that 
the  concentration  of  each  subset  and,  in  turn,  the  ratio  between  subsets  can  be  predicted. 

7.3  Predicting  Combined  Injury  in  an  Urban  Environment 

To  improve  operational  utility  of  the  hematopoietic  models,  the  suite  of  models  will  even¬ 
tually  be  linked  to  other  tools  which  provide  environment  inputs  for  modern  scenarios  that 
include  urban  environments.  Current  tools,  such  as  HPAC  and  NucFast  calculate  free-in-air 
radiation  dose  as  well  as  the  incident  thermal  fluence  at  specihc  locations  in  urban  environ¬ 
ments.  Future  work  will  involve  correlating  the  direct  thermal  fluence  from  the  blast  and 
secondary  hres  with  %TBSA  burn.  BURNSIM,  a  tool  which  models  thermal  injury  based  on 
the  heat  flux,  may  be  useful  in  this  task  (Knox  2007).  This  tool  provides  information  on  the 
burn-depth  using  heat  flux  and  various  types  of  clothing  as  inputs.  By  setting  a  threshold 
burn-depth  level  and  taking  into  account  the  fraction  of  exposed/covered  skin,  it  may  be 
possible  to  use  this  tool  to  relate  thermal  fluence  to  %TBSA  burned.  With  the  integration 
of  these  models  with  geospatial  population  data,  more  detailed  and  precise  combined  injury 
calculations  in  urban  scenarios  will  be  possible. 
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8  Conclusion 


Models  of  thrombopoiesis,  granulopoiesis,  and  lymphopoiesis  have  been  developed  to  simu¬ 
late  the  effects  of  burn  injury  on  hematopoietic  cell  kinetics  and  integrated  with  radiation 
effects  models  to  make  combined  injury  predictions.  The  model  uses  acute  radiation  expo¬ 
sure  and  percent  surface  area  burned  as  inputs  and  provides  predictions  of  blood  cell  and 
progenitor  cell  counts  over  time.  Each  model  has  been  optimized  and  validated  against 
human  data  for  burn  and  radiation  exposure  separately.  Model  outputs  were  compared  to 
combined  injury  data  from  animals  to  verify  that  trends  were  accurately  predicted.  Further 
validation  will  be  required  if  new  data  from  patients  with  combined  injury  becomes  available. 
Outputs  from  these  models  (e.g.,  length  of  cytopenia  or  depth  of  nadirs)  will  be  correlated 
with  mortality  to  establish  a  predictive  model  of  mortality  risk.  This  information  will  assist 
in  disaster  preparedness  and  medical  resource  planning. 
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Appendix  A  Equations  for  Hematopoietic  Models 

A.l  Equations  for  Thrombopoiesis  Model 

Below  is  the  complete  list  of  equations  for  the  human  model  of  thrombopoiesis. 


Cell  concentrations  after  acute  radiation  exposure  (dose  D)  at  time  t  =  0: 
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A. 2  Equations  for  Granulopoiesis  Model 

Below  is  the  complete  list  of  equations  for  the  human  model  of  granulopoiesis. 

Cell  concentrations  after  acnte  radiation  exposure  (dose  D)  at  time  f  =  0: 
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Change  in  cell  concentrations  with  time: 
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S  is  the  fraction  of  skin  surface  area  affected  by  burn  (0  <  S'  <  1),  and 
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A. 3  Equations  for  Lymphopoiesis  Modei 

Below  is  the  complete  list  of  equations  for  the  human  model  of  lymphopoiesis. 

Cell  concentrations  after  acnte  radiation  exposure  (dose  D)  at  time  t  =  0: 


xf  =  xf  (l  -  (1  -  )"*)  {i  =  1,  2,  3)  (62) 

xf  =  x“"(l-e-^/^°)"‘  (*  =  1,2,3)  (63) 

Change  in  cell  concentrations  with  time: 

xf  =  Bxf  -  7xf  (64) 

xf  =  7xf  -  bxf  (65) 

xf  =  bxf  -  %l)xf  (66) 

xf  =  -^xf  (i  =  l,2,3)  (67) 

where 


1  +  /3{xi  +  xf  +  92{x2  +  xf  +  Ofxs  +  xf)) 

^  =  ^o(l  +  &3^"^e-“^*) 


and  S  is  the  fraction  of  skin  surface  area  affected  by  burn  (0  <  S'  <  1). 


(68) 

(69) 
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Appendix  B  Murine  Thrombopoiesis  and  Granulopoiesis 

Models 

In  this  section  the  development  of  the  murine  thrombopoiesis  and  granulopoiesis  models 
is  discussed.  The  equations  for  the  models  are  the  same  as  those  given  in  Appendix  A. 
However,  the  parameter  values  in  the  models  were  changed  using  data  from  rodent  studies. 

B.l  Methods 

Parameter  values  were  optimized  using  differential  evolution  optimization  with  the  DEoptim 
library  in  R  (R  Core  Team  2015;  Mullen  et  ah  2011;  Ardia  et  ah  2015).  To  run  the  opti¬ 
mization,  an  upper  and  lower  bound  for  each  parameter  to  be  optimized  were  specihed  along 
with  a  cost  value.  For  our  purposes,  we  set  this  cost  value  equal  to  the  sum  of  the  squared 
residuals.  The  optimization  was  run  using  normalized  versions  of  the  data  (i.e.,  data  points 
were  divided  by  a  constant  baseline  value).  For  doses  given  in  roentgen,  a  0.0096  Gy  per 
roentgen  conversion  factor  was  used.  The  optimized  parameters  and  the  bounds  imposed  are 
given  in  Table  B.l  and  B.2  for  the  thrombopoiesis  and  granulopoiesis  models,  respectively. 
The  models  were  first  optimized  to  murine  radiation  data.  Next,  burn-dependent  parameters 
were  determined  by  optimizing  the  models  to  murine  burn  data. 

Table  B.l:  Parameter  bounds  used  during  the  murine  thrombopoiesis  model 
optimizations. 


Parameter  Lower  Bound  Upper  Bound 


7  (d-^) 

0.5 

1.9 

92 

0.001 

1 

A 

0 

6 

(Gy) 

0.5 

3 

D'i  (Gy) 

0 

30 

bo 

0 

20 

ao  (d“^) 

0.001 

10 

Co 

0.1 

20 

bf 

0 

20 

a/  (d  ^) 

0 

10 

c/ 

0.1 

20 
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Table  B.2:  Parameter  bounds  used  during  the  murine  granulopoiesis  model  op¬ 
timizations. 


Parameter  Lower  Bound  Upper  Bound 


7]  (d  ^) 

0 

2 

ip  (d-i) 

0 

15 

m 

0 

1000 

(Gy) 

0 

0.899 

Dr  (Gy) 

0.9 

10 

^2°  (Gy) 

0 

15 

bo 

0 

10 

ao  (d"^) 

0.0001 

0.1 

b2 

5 

100 

02  (d"^) 

0.05 

10 

bs 

0 

20 

03  (d“^) 

0 

10 

A 

0 

1000 

k 

0 

10 

B.2  Thrombopoiesis  Murine  Model  Results 

Table  B.3  provides  the  finalized  parameter  values  for  the  murine  thrombopoiesis  model.  Of 
the  optimized  parameter  values,  6*2  approached  the  lower  bound  of  0.001,  A  approached 
the  upper  bound  of  6.0,  and  bf  approached  the  upper  bound  of  20.  No  other  optimized 
parameters  approached  the  bounds  specified.  Sources  for  parameters  that  were  not  optimized 
are  given  in  Table  B.3.  The  derivations  of  many  of  these  parameter  values  are  described 
elsewhere  (Wentz  2014a;  Wentz  et  al.  2014a),  while  the  derivations  of  (5o  and  '0  for  the 
murine  thrombopoiesis  model  are  given  below. 

There  are  three  stages  of  megakaryocyte  maturation:  megakaryoblast,  promegakaryocyte, 
and  granular  megakaryocyte.  The  maturation  time  through  all  three  compartments  in  mice 
is  approximately  68  h  (Ebbe  1971).  However,  rodent  platelet  counts  following  radiation 
exposure  at  low  doses  typically  do  not  decline  until  4  days  post  irradiation  (Ebbe  et  al. 
1970).  This  suggests  that  there  is  a  fourth  sub-compartment  in  the  post-mitotic  compartment 
which  represents  non-mitotic  megakaryocyte  progenitors.  Using  4  days  as  the  maturation 
time  leads  to  a  normal  maturation  rate,  Sq,  of  0.25  d“^. 

Platelet  lifespans  vary  between  mice  strain;  one  experiment  noted  the  platelets  in  C3H 
mice  have  a  shorter  lifespan  (4.8  d)  than  platelets  in  C57BL  mice  (5.7  d)  (Manning  et  al. 
1997).  Other  experiments  reported  platelet  lifespans  between  3.1  and  3.3  d  (Baker  et  al. 
1997).  In  the  model  here,  the  expected  platelet  lifespan  is  4  d,  a  generally  accepted  value 
(Wickrema  et  al.  2009).  Thus,  ■0  is  1/4  or  0.25  d“h 
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Table  B.3:  Biological  parameter  values  for  murine  thrombopoiesis  model. 


Parameter 

Description 

Value 

Section  (References) 

a 

Maximum  rate  of  mitotic  progenitor  cell  repop¬ 
ulation 

2.0  d-i 

Andreeff  et  al.  2000* 

7 

Rate  of  mitotic  progenitor  cell  maturation 

1.60  d-i 

Optimized 

<5o 

Rate  of  MK  maturation 

0.25  d-i 

Ebbe  1971 

V' 

Rate  of  platelet  decay 

0.25  d-i 

Wickrema  et  al.  2009 

n 

Number  of  MK  subcompartments 

10 

Same  as  human  model* 

m 

Number  of  platelet  subcompartments 

9 

Same  as  human  model* 

02 

Decay  rate  of  mediator  due  to  X2  cells  relative 
to  decay  rate  due  to  xi  cells 

0.001 

Optimized 

03 

Decay  rate  of  mediator  due  to  X3  cells  relative 
by  decay  rate  due  to  xi  cells 

0.03  02 

Kuwaki  et  al.  1998; 
Broudy  et  al.  1997* 

^min 

Minimum  rate  of  MK  maturation 

(4/3)S 

Same  as  human  model* 

^max 

Maximum  rate  of  MK  maturation 

4<5 

Same  as  human  model* 

A 

Strength  of  X3  feedback  on  immature  MK  mat¬ 
uration 

6.0 

Optimized 

a 

Number  of  platelets  per  MK 

3000 

Smirnova  2010; 

Deutsch  et  al.  2013 

D°i 

Determines  fraction  of  xi  targets  that  are  hit  by 
radiation  at  dose  D 

1.37  Gy 

Optimized 

ni 

Number  of  targets  per  xi  cell 

1 

Same  as  human  model* 

D°2 

Determines  fraction  of  X2  targets  that  are  hit  by 
radiation  at  dose  D 

4.65  Gy 

Optimized 

U2 

Number  of  targets  per  X2  cell 

4 

Same  as  human  model* 

Rate  of  moderately  damaged  cell  death 

1.0  d-i 

Smirnova  2010* 

bo 

Determines  max  burn  effect  on  repopulation 
rate 

5.13 

Optimized 

ao 

Determines  duration  of  burn  effect  on  repopula¬ 
tion  rate 

0.043  d-i 

Optimized 

Co 

Determines  how  run  size  affects  repopulation 
rate  change 

1.3 

Optimized 

bf 

Determines  max  burn  effect  on  MK  maturation 

20 

Optimized 

a/ 

Determines  duration  of  burn  effect  on  MK  mat¬ 
uration 

0.00017 

Optimized 

c/ 

Determines  how  burn  size  affects  MK  matura¬ 
tion  rate  change 

1.6 

Optimized 

do 

Determines  max  effect  of  burn  on  platelet  decay 

0.476 

Wentz  2014a 

di 

Determines  duration  of  effect  of  burn  on  platelet 
decay 

0.143  d-i 

Wentz  2014a 

*See  Wentz  et  al.  2014a  for  description  of  how  parameter  value  was  derived. 


Radiation  Optimization  Results 

Figure  B.l  shows  comparisons  of  the  murine  thrombopoiesis  model  with  data  used  for  opti¬ 
mization  of  the  baseline  and  radiation  portion  of  the  model. 
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Platelet  concentration  (normalized)  Platelet  concentration  (normalized) 


(a)  1.92  Gy  (Ebbe  et  al.  1970)  (b)  2.88  Gy  (Ebbe  et  al.  1970) 


(c)  6  Gy  (Ghosh  et  al.  2012)  (d)  7  Gy  (Satyamitra  et  al.  2012; 

Neumedicines,  Inc.) 

Figure  B.l:  Murine  thrombopoiesis  model  compared  with  radiation  optimization 
data  (mean  ±  standard  error). 
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Radiation  Validation  Results 

Figure  B.2  show  comparisons  of  the  murine  thrombopoiesis  model  with  validation  data. 


(a)  0.96  Gy  (Cohn  et  al.  1956) 


(b)  1.92  Gy  (Ebbe  et  al.  1970) 


(c)  2.4  Gy  (Cohn  et  al.  1956)  (d)  2.88  Gy  (Ebbe  et  al.  1970) 

Figure  B.2:  (1  of  4)  Murine  thrombopoiesis  model  compared  with  validation 
data  (mean). 
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Platelet  concentration  (normalized)  Platelet  concentration  (normalized) 


(e)  3  Gy  (Singh  et  al.  2009)  (f)  3.84  Gy  (Gohn  et  al.  1956;  Ebbe  et  al. 

1970) 


(g)  4.8  Gy  (Gohn  et  al.  1956)  (h)  6  Gy  (Faltynek  et  al.  1992) 

Figure  B.2:  (2  of  4)  Murine  thrombopoiesis  model  compared  with  validation 
data  (mean). 
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(i)  6  Gy  (Cohn  et  al.  1956) 


(k)  7  Gy  (Satyamitra  et  al.  2011;  Singh 
et  al.  2009) 

Figure  B.2: 
data  (mean). 


(j)  6.72  Gy  (Cohn  et  al.  1956) 


(1)  7.5  Gy  (Davis  et  al.  2010) 


(3  of  4)  Murine  thrombopoiesis  model  compared  with  validation 
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(m)  7.68  Gy  (Cohn  et  al.  1956) 

Figure  B.2:  (4  of  4)  Murine  thrombopoiesis  model  compared  with  validation 
data  (mean). 


Burn  Optimization  Results 

Figure  B.3  shows  comparisons  of  the  thrombopoiesis  model  with  data  used  to  optimize  the 
burn-dependent  parameters  in  the  model. 
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(a)  15%  TBSA  (Davis  et  al.  1955) 


(b)  20%  TBSA  (Wallner  et  al.  1984) 


Time  following  burn  (days) 


(c)  20%  TBSA  (Kalmaz  et  al.  1991) 


(d)  25%  TBSA  (Fujimi  et  al.  2006) 


Figure  B.3:  (1  of  2)  Murine  thrombopoiesis  model  compared  with  burn  opti¬ 
mization  data  (mean  ±  standard  error). 
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Platelet  concentration  (normalized) 


(e)  25%  TBSA  (Davis  et  at.  1955) 


(f)  30%  TBSA  (Eurenius  et  al.  1972) 


Figure  B.3:  (2  of  2) 
mization  data  (mean 


Murine  thrombopoiesis  model  compared  with 
±  standard  error). 


burn  opti- 


B.3  Granulopoiesis  Murine  Model  Results 

Table  B.4  provides  the  finalized  parameter  values  for  the  murine  granulopoiesis  model.  Of 
the  optimized  parameter  values,  -0  approached  the  upper  bound  of  15  d“^,  oq  approached 
the  lower  bound  of  0.0001  d“^,  and  02  approached  the  lower  bound  of  0.05  d“^.  All  other 
parameter  values  were  optimized  to  a  value  within  the  bounds  specihed.  Sources  for  pa¬ 
rameter  values  that  were  not  optimized  are  given  in  Table  B.4.  The  derivations  of  many  of 
these  parameter  values  are  described  elsewhere  (Wentz  2014b;  Wentz  et  ah  2014a),  while  the 
derivations  of  7,  k,  6,  I,  r,  and  v  for  the  murine  granulopoiesis  model  are  model  are  given 
below. 

Compartment  transition  rates  were  estimated  using  experimental  data  in  mice.  The 
values  for  7  and  k  were  taken  directly  from  the  literature  (Pillay  et  ah  2010).  The  transit  time 
through  the  mitotic  compartment  was  0.7  d,  while  the  circulating  lifespan  of  the  neutrophils 
was  0.75  d.  This  leads  to  rates  of  1.43  d“^  and  1.33  d“^  for  7  and  n,  respectively.  The  transit 
time  through  the  post- mitotic  compartment  was  measured  as  1.6  d  (Pillay  et  ah  2010).  In 
the  model,  this  rate  is  determined  by  6,  m,  and  I  where  the  rate  of  transition  F  is 


^  1  +  mxl/xl 
1  + Ixl/xl 


(70) 


Where  6  is  the  maximum  rate  of  transition,  and  5m /I  is  the  minimum.  Assuming  the  data 
was  obtained  at  equilibrium  (X3  =  X3)  and  Fg  is  the  transition  rate  at  equilibrium: 


F,  =  5 


1  +  m 

l  +  l 


(71) 


Since  the  post-mitotic  transition  time  was  measured  as  1.6  d,  F^  =  0.625  d  ^  and 

5  =  0.625^^  (72) 

1  -|-  m 

A  study  in  mice  examined  how  injected  G-CSF  altered  the  post-mitotic  transition  time. 
The  transition  time  was  found  to  be  reduced  by  55%  (Uchida  et  ah  1992).  Assuming  this 
represents  the  minimum  transition  time  or  maximum  rate  (5),  we  can  derive  a  relationship 
between  m  and  h 


1  1 

(73) 

(74) 

^  1  +  /  =  2.22(1  +  m) 

(75) 

^  1  =  1.22  + 2.22m 

(76) 

The  time  of  abortive  elevation  is  determined  using  the  following  equation: 

Atae  =  T  -vD  (77) 
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where  D  is  the  radiation  dose,  r  and  v  were  determined  throngh  examination  of  the  time 
of  abortive  rise  in  experimental  data.  At  D=2  Gy,  the  abortive  rise  occnrred  approximately 
on  day  10  {Atae  =  10)  (Romero- Weaver  et  ah  2012;  Romero- Weaver  et  ah  2013).  At  D=7 
Gy,  the  abortive  rise  occnrred  approximately  on  day  7  (Nenmedicines,  Inc.).  Using  these 
two  time  points,  valnes  of  11.2  for  r  and  0.6  for  v  were  determined. 
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Table  B.4:  Biological  parameter  values  for  murine  granulopoiesis  model. 


Parameter 

Description 

Value 

Source/Reference 

a 

Maximum  rate  of  mitotic  progenitor  cell 
repopulation 

2.0  d-i 

Andreeff  et  al.  2000* 

7 

Rate  of  mitotic  progenitor  cell  maturation 

1.43  d-i 

Pillay  et  al.  2010 

<5 

Determines  post-mitotic  maturation  rate 

0.63  ^  d-i 

Pillay  et  al.  2010 

m 

Determines  post-mitotic  maturation  rate 

11.54 

Optimized 

1 

Determines  post-mitotic  maturation  rate 

1.22  -h  2.22m 

Uchida  et  al.  1992 

K 

Rate  of  transition  from  blood  to  tissue 

1.33  d-i 

Pillay  et  al.  2010 

V' 

Rate  of  granulocyte  decay  from  tissue 

15.0  d-i 

Optimized 

V 

Rate  of  weakly  damaged  cell  death  (when 
t  >  Atae) 

0.28  d-i 

Optimized 

02 

Decay  rate  of  mediator  due  to  X2  cells  relative  to 
decay  rate  due  to  xi  cells 

2.24 

Shinjo  et  al.  1995* 

03 

Decay  rate  of  mediator  due  to  X3  cells  relative  to 
decay  rate  due  to  xi  cells 

8.45 

Shinjo  et  al.  1995* 

04 

Decay  rate  of  mediator  due  to  X4  cells  relative  to 
decay  rate  due  to  xi  cells 

03 

Assumed 

Dl 

Determines  fraction  of  damaged  xi  cells 

0.67  Gy 

Optimized 

DOOo 

Determines  ratio  of  weakly  damaged  to  damaged 
xi  cells 

1.14  Gy 

Optimized 

Determines  fraction  of  damaged  X2  cells 

8.12  Gy 

Optimized 

Dl 

Determines  fraction  of  damaged  X3  cells 

Dl 

Assumed 

Dl 

Determines  fraction  of  damaged  X4  cells 

Dl 

Assumed 

T 

Determines  time  of  abortive  rise: 

Atae  =  T  -vD 

11.2  d 

See  text 

V 

Determines  time  of  abortive  rise: 

Atae  =  T  -vD 

0.6  d  Gy-i 

See  text 

rii 

Number  of  targets  per  cell  Xi 
(z  =  l,2,3,4) 

1 

Assumed 

M 

Rate  of  damaged  cell  death 

1.0  d-i 

Smirnova  2010* 

bo 

Determines  max  burn  effect  on  repopulation  rate 

0.381 

Optimized 

ao 

Determines  duration  of  burn  effect  on  repopulation 
rate 

0.0001  d-i 

Optimized 

^2 

Max  relative  burn  effect  on  bone  marrow  release 

6.373 

Optimized 

02 

Determines  duration  of  burn  effect  on  bone 
marrow  release 

0.050  d-i 

Optimized 

03 

Max  relative  burn  effect  on  transition  to  tissue 

0.055 

Optimized 

03 

Determines  duration  of  burn  effect  on  transition  to 
tissue 

0.036  d-i 

Optimized 

br 

Determines  max  burn  effect  on  marginated  ratio 

3 

Wentz  2014b 

ar 

Determines  duration  of  burn  effect  on  marginated 
ratio 

0.106  d-i 

Wentz  2014b 

A 

Determines  how  changing  burn  size  affects  the 
system 

962 

Optimized 

k 

Determines  how  changing  burn  size  affects  the 
system 

0.336 

Optimized 

*See  Wentz  et  al.  2014a  for  description  of  how  parameter  value  was  derived  from  source. 
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Radiation  Optimization  Results 

Figure  B.4  shows  comparisons  of  the  murine  granulopoiesis  model  with  data  used  for  opti¬ 
mization  of  the  baseline  and  radiation  portion  of  the  model. 


(a)  2  Gy  (Diamonds:  Romero- Weaver 
et  ah  2012;  Squares:  Romero- Weaver 
et  ah  2013) 


(b)  3.84  Gy  (Brecher  et  al.  1948) 


Figure  B.4:  (1  of  2)  Murine  granulopoiesis  model  compared  with  radiation  op¬ 
timization  data  (mean  ±  standard  error). 
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(c)  5  Gy  (Palmer  et  al.  2011) 


(d)  6  Gy  (Ghosh  et  al.  2012) 


(e)  7  Gy  (Diamonds:  Satyamitra  et  al. 

2012;  Squares:  Neumedicines,  Inc.) 

Figure  B.4:  (2  of  2)  Murine  granulopoiesis  model  compared  with  optimization 
data  (mean  ±  standard  error). 
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Radiation  Validation  Results 

Figure  B.5  show  comparisons  of  the  murine  granulopoiesis  model  with  validation  data. 


(b)  4.56  Gy  (Smith  et  al.  1954) 


10 


20 


30 


Time  following  radiation  (days) 

(c)  5.76  Gy  (Smith  et  al.  1954) 


(d)  7  Gy  (Diamonds:  Singh  et  al.  2009; 
Squares:  Satyamitra  et  al.  2011) 


Figure  B.5:  (1  of  2)  Murine  granulopoiesis  model  compared  with  validation  data 
(mean). 
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Granulocyte  concentration  (normalized) 


(e)  7.5  Gy  (Davis  et  al.  2010)  (f)  8  Gy  (Kulkarni  et  al.  2013) 


Radiation  dose  (Gy) 

(g)  7  days  post  exposure  (Tan  et  al.  2006; 
errors  bars  represent  standard  error) 

Figure  B.5:  (2  of  2)  Murine  granulopoiesis  model  compared  with  validation  data 
(mean). 
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Burn  Optimization  Results 

Figure  B.6  shows  comparisons  of  the  granulopoiesis  model  with  data  used  to  optimize  the 
burn-dependent  parameters  in  the  model. 


(a)  10%  TBSA  (Gruber  et  ah  1989)  (b)  15%  TBSA  (Palmer  et  ah  2011) 

Figure  B.6:  (1  of  2)  Murine  granulopoiesis  model  compared  with  burn  optimiza¬ 
tion  data  (mean  ±  standard  error). 
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(c)  20%  TBSA  (Gruber  et  al.  1989)  (d)  20%  TBSA  (Wallner  et  al.  1984) 


(e)  30%  TBSA  (Gruber  et  al.  1989) 


(f)  60%  TBSA  (McManus  1983) 


Figure  B.6:  (2  of  2)  Murine  granulopoiesis  model  compared  with  burn  optimiza¬ 
tion  data  (mean  ±  standard  error). 
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Abbreviations,  Acronyms,  and  Symbols 


ARA 

d 

DTRA 

FME 

G-CSF 

Gy 

h 

HENRE 

HPAC 

IL-1 

IL-6 

MCMC 

MK 

NA 

/iL 

NK 

ODE 

R 

S 

SSR 

TBSA 

TPO 


Applied  Research  Associates,  Inc. 
day 

Defense  Threat  Reduction  Agency 

flexible  modeling  environment  for  modeling,  sensitivity,  and  Monte  Carlo 
analysis  (R  plug-in) 

Granulocyte  colony  stimulating  factor 

gray 

hour 

Health  Effects  from  Nuclear  and  Radiation  Environments  (HENRE) 
Hazard  Prediction  and  Assessment  Capability 
interleukin  1 
interleukin  6 

Markov  Chain  Monte  Carlo 

megakaryocyte 

not  available 

microliter 

natural  killer 

ordinary  differential  equation 

software  programming  language  for  statistical  computing 

surface  area  burn  size  fraction 

sum  of  squared  residuals 

total  body  surface  area 

thrombopoietin 
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